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Preface 



This thesis is submitted in partial fulfillment of the PhD degree at the Niels 
Bohr Insitute in Denmark. The thesis concludes three years of work, of 
which two years were dedicated to scientific work and one year to lecturing 
and various other duties. The PhD project was carried out under supervi- 
sion of Professor Ake Nordlund, in the Computational Astrophysics group 
at the Niels Bohr institute in Copenhagen, Denmark. As part of the PhD 
education, I spent 4 months at the National Space Science and Technology 
Center (NSSTC) in Huntsville Alabama, USA. Here I collaborated with the 
gamma-ray group under Dr. Gerald J. Fishman and in particular with Dr. 
Ken-Ichi Nishikawa. 

The backbone of the thesis consists of three papers published in Ap JL (the 
Astrophysical Journal Letters), plus recent work on large-scale 2D plasma 
simulations and generation of synthetic spectra that has not yet been sub- 
mitted as scientific papers. The three papers are included as chapters here, 
with minor extensions compared to the published versions. Co-author state- 
ments are attached at the end of the thesis. 

The aim has been to write a thesis that both students and experienced 
scientist may benefit from reading. This is a difficult undertaking and thus 
you may find some parts either too trivial or too complex. In the latter case, 
do not worry: Collisionless plasma shocks are indeed complicated non-linear 
systems, and really a contradiction in terms. 

I hope that you will enjoy and appreciate the work I present in this thesis. 



Christian Hededal 

Copenhagen, May 2005 



The illustration on the cover is a ray-traced image of data, taken directly from our 
simulations. It shows plasma filaments that are generated by the Weibel two-stream 
instability, as a relativistic plasma is travelling directly towards the observer. 



Abstract 



In this thesis I present the results of three-dimensional, self-consistent particle- 
in-cell simulations of collisionless shocks. 

The radiation from afterglows of gamma-ray bursts is generated in colli- 
sionless plasma shocks between a relativistic outflow and a quiescent circum- 
burst medium. The two main ingredients responsible for the radiation are 
high-energy, non-thermal electrons {N{'-f)d'~f oc 7~^) and a strong magnetic 
field. Fermi acceleration is normally believed to be responsible for the accel- 
eration of the electrons. Fermi acceleration has been employed extensively in 
Monte Carlo simulations, where it operates in conjunction with certain as- 
sumptions about the scattering of particles and the structure of the magnetic 
field. The mechanism has, however, not been conclusively demonstrated to 
occur in ab initio particle simulations and also faces additional problems. 
Furthermore, the requirement of a strong magnetic field in the shock region 
indicates that the magnetic field is generated in situ in. 

In this thesis, I argue that in order to make the right conclusions about 
gamma-ray burst and afterglow parameters, it is crucial to have a firm under- 
standing of collisionless shocks: How are electrons accelerated, what is the 
topology of the generated magnetic field, and how do these two aspects af- 
fect the resulting radiation. Thus, the main goal of the work I present in this 
thesis has been to expand our knowledge about the microphysics of collision- 
less plasma shocks. To accomplish this, a self-consistent, three-dimensional 
particle-in-cell computational code has been utilized. The simulation tool 
works from first principles by solving Maxwell's equations for the electromag- 
netic field, consistently coupled to the momentum equation for the charged 
particles. 

In the experiments, I study the collision of two plasma populations trav- 
elling at relativistic velocities. When the plasma populations are initially 
unmagnetized or weakly magnetized, the Weibel two-stream instability gen- 
erates a magnetic field in the shock ramp with strengths up to percents of 
equipartition with the plasma ions. The nature of the magnetic field is pre- 
dominantly transverse to the plasma flow. The transverse coalescence scale 



iii 



is comparable to the ion skin depth whereas the parallel scale extends up to 
hundreds of ion skin depths. A spatial Fourier decomposition of the magnetic 
field shows that the structures follow a power-law distribution with negative 
slope. 

The experiments also reveal a new non-thermal electron acceleration 
mechanism, which differs substantially from Fermi acceleration. Accelera- 
tion of electrons is directly related to the formation of ion current channels 
by the non-linear Weibel two-stream instability. This links particle acceler- 
ation closely together with magnetic field generation in coUisionless shocks. 
The resulting electron spectrum consists of a thermal component and a non- 
thermal component at high energies. In an experiment with a bulk Lorentz 
factor of r = 15, the non-thermal tail has the power-law index p = 2.7. 

Finally, I have developed a tool that generates synthetic radiation spectra 
from the experiments. The radiation is calculated directly, by tracing a large 
number of electrons in the generated magnetic field, and thus continuous the 
line of work from first principles. Numerous tests show that the radiation 
tool successfully reproduces synchrotron, bremsstrahlung and undulator ra- 
diation from small-angle deflections. I then go on to perform a parameter 
study of three-dimensional jitter radiation. Using the tool on particle-in-cell 
experiments of coUisionless shocks I find that the radiation spectrum from 
particles in a randomized magnetic field is not fully consistent with radiation 
from particles in shock-generated magnetic field, even when the two have the 
same statistical properties. 

In experiments where magnetic field generation and particle acceleration 
arise as natural consequences of the Weibel two-stream instability, the result- 
ing radiation spectrum is consistent with observations. In simulations of a 
coUisionless shock that propagates with bulk Lorentz factor F = 15 through 
the interstellar medium, I find that the radiation spectrum peaks around 
lO^^Hz. Above this frequency, the spectrum follows a power-low F oc 
with /3 = 0.7. Below the peak frequency, the spectrum follows a power law 
F oc z/" with a ~ 2/3. This is steeper than the standard synchrotron value 
of 1/3 and more compatible with observations. 

I conclude that strong magnetic field generation (e^ ~ 0.01 — 0.1), non- 
thermal particle acceleration, and the emission of radiation with properties 
that are consistent with GRB afterglow observations are all unavoidable con- 
sequences of the collision between two relativistic plasma shells. 
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Units and Conventions 



The field of gamma-ray bursts is multi-disciplinary in the sense that it covers 
the electromagnetic spectrum from radio to gamma-ray, length scales from 
microphysical to cosmological and velocities from newtonian to highly rela- 
tivistic. It thus connects many branches of physics and astrophysics, with 
all their traditions with regard to formalisms, phenomenologies and units. 
Some people are almost religious about units. I'm not, as long as the use of 
units is consistent. I have chosen to use SI- units throughout the thesis. To 
those readers that are mostly familiar with the use of Gaussian cgs units in 
electrodynamics, the following conversion table may be of help: 

SI cgs 
Replace by 

Aio 47r/c^ 
B B/c 

Throughout the thesis I use 7 for Lorentz of individual particles and F for 
bulk flows, e.g. jet Lorentz factor. On many occasions, I use the terms "par- 
allel" and "perpendicular". If nothing else is stated, this refers a direction 
relative to propagation direction of the shock. 

Since the units in PIC codes are often re-scaled (and our code is no ex- 
ception) it is normal to measure time and length in terms of the typical 
plasma parameters that govern the physical processes in a plasma. Time 
is often given in units of one over the electron plasma frequency Upe = 
{n(,q^ / (mgeo))^^^, where rig is the electron plasma density, q is the unit charge, 
TTLf. is the electron mass, and eo is the electric vacuum permittivity. Lengths 
are typically given in units of skin depths 5 = c/upe where c is the speed 
of light in vacuum. If not otherwise stated, the units in the figures should 
be taken as arbitrary units. However, on some occasions I make an effort 
to scale the results into real-space values. In this case, the units are clearly 
marked on the axis. 
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Chapter 1 
Introduction 



1.1 The "early" history of gamma-ray bursts 

Gamma-ray bursts are the largest explosions we know of in the 
universe after the Big Bang 

If you have seen one gamma-ray burst... you have seen one 
gamma-ray burst! 



These two statements are among the most repeated in the history of gamma- 
ray bursts. They very nicely cover what the fuss about gamma-ray bursts 
is all about. With their extreme brightness, gamma-ray bursts have the 
potential of being used as lighthouses, shining from the far and dark ages 
of the universe. At the same time, gamma-ray bursts apparently come in a 
large number of colors and flavors, and thus their origin is still a puzzle after 
many years in scientific focus. 

The existence of gamma-ray bursts (GRBs) came into human knowledge 
at the end of the 1960s. As a product of the nuclear arms race, the USA 
had launched a series of gamma-ray nuclear blast wave detectors — the 
VELA satellites. The aim of the VELA satellites was to make sure that the 
USSR did not break the nuclear test ban treaties with secret nuclear tests in 



the up per atmosphere and in space. Testing the satellites, iKlebesadel et al. 



( 19731 ) found unidentified spikes in the data. It was easily realized that the 
signals were not from nuclear tests. Using the timing offset from several 
satellites, it was possible to make a crude triangulation and place the origin 
of the gamma-rays outside our Solar system. Distributed randomly in the 
sky, the positions indicated that the bursts were either from an extended 
galactic halo, or were an extra-galactic phenomenon. 
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In 1991, the Compton Gamma Ray Observatory (CGRO) was launched, 
carrying the Burst and Transient Source Experiment (BATSE). Severa l thou- 



sand detections during the 1990s, isotropically distributed over the sky (jMeegan et al 
1993), still lef t two possibilities for the origin of the GRBs. Either they were 



cosmological (jPaczvhski. Il995t ) or th ey were from a very extended spherical 
halo of the Milky way (lLambl . ll995h . Which of the two was determined in 
1997, following the launch of the BeppoSAX satellite. BeppoSAX was able 
to rapidly locate the position of GRB 970228 (970228 for 1997, February 28). 
This triggered a multiwavele ngth campaign resulting in the detection of an 
x-ray afterglow (Costa, 1997), and an optical afterglow ( van Paradiid . 19971 ) 
within the error-box. Using the Hubble Space Telescope, the or i gin of the 
burst was found to lie in a galaxy at cosmological distance (jSahd . Il997l ). 
The sp ectrum of the host g alaxy to GRB 970228 showed prominent emis- 



sion lines (jBloom et al.L 1200 1[ ). From these, the redshift of the galaxy was 
found to be 2; = 0.695. With the given flux and enormous distance, the 
energy of the GRB was estimated to be as large as 10^^ J (10^^ erg), many 
thousand times stronger than any previous known type of astrophysical ex- 
plosion. Moreover, the duration over which the energy was released was 
apparently only a few seconds. The great variability in the burst suggested 
that this huge amount of e nergy was released within a volume a few 1000 
km in radius. Ruderman ( 1975| ) had realized that such a scenario would 
inevitably lead to a compactness problem. The fireball would be extremely 
optically thick with respect to pair production and this would not allow us to 
observe the high-energy, non-thermal photon tail. This problem was solved 
by suggesting tha t the emitting surface was e j ected with a highly relativis- 
tic bulk velocity ( Goodman . 19861 : Paczvnski . 19861 : Piran, 1996[ ). Fireball 
expansion speeds comparable to the speed of light where later o bservation- 
ally co nfirr ned from changes in ra dio scintillation of GRB 970508. [Goodman 
(1997) and Waxman et al. ( 1998[ ) estimated the size of the fireball to lO^^m, 
only four weeks after the trigger. 



In 1998, supernova ib)b)»sw was louna w itnm tne error box or (jtin y»U4:^£) 
( Galama et al.1 . 19981 : Kulkarni et al. . 19981 ) and in 2003, the Supernova-GRB 
connection was unambiguous established by the discovery of a clear super- 
nova lig ht-curve bump and spectral signature in the optical afterglow of GRB 



jernova 1998sw was found within the error box of GRB 980425 



030329 (iHiorth et al.1 



fana s 
200I 



Stanek et al.1 12003> ) (see Fig. EH). 



1.2 The general picture 

After roughly 40 years with GRBs in the scientific spotlight, we are converg- 
ing towards a working theory for GRBs. 
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Figure 1.1: Left: Spectral evolution of the combined optical flux density of 
GRB 030329, the associated SN2003dh, and its host galaxy {colored lines) com- 
pared to the spectrum of SN1998bw (dashed line). Right: R-band light curve for 
GRB041006. A typical GRB power-law decay of the optical afterglow is shown 
as the dashed line and the SN1998bw light curve extended by a fac tor of 1.38 is 



shown as the dot ted curve. Together they fit the observations. From lHiorth et al 



(|2003l ^ {left) and lStanek et a].r(|200,'j ^ {right). 



There appears to be a general consensus in the scientific community about 
the fireball internal-external shock model, in which the gamma-ray burst and 
the subsequent afterglow radiation is created by dissipation of collisionless 
plasma shocks. Independent of the true nature of the progenitor, the ex- 
tremely large amount of energy deposited in a very small volume inevitably 
creates a highly energetic outflo w that will interact with the surrounding 
medium ( Shemi and Piran . 199fl[ ). 

For a detailed discussion, several good review papers exist fe.g., Fishman and Meeganlll995L 



- ,.pers exist C . 

PiraJl999[lva,n Paradiis et a llEiil iMas/irolEiml iMeszarojEnilZhang et abkonT 
and lPiranllinnibh . 



To explain the origin and variability of the prompt gamma emission, it 
was suggested that the progenitor may expel multiple plasma shells with 
differe nt energies. The shells h e at up in shocks when th ey overtake each 
other ( Rees and Meszarosl . 1994 : Paczvnski and Xu . 1994 ). At later times, 
an external shock forms as the ejecta blasts through the external medium. 
The external medium can either be the interstellar medium or a progenitor 
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wind. This shock heats up the external plasma an d creates the afterglow 
( Rees and Meszarosl . 19921 : Meszaros and Reei . 19931 ). 

As the fireball expands into the external medium it sweeps up external 
matter and is e xpected to eventually approa ch the self similar Blandford- 
McKee solution f)Blandford and McKeel . Il97fit) 



This is a blast wave solution 



analogues to the non-rela tivistic Sedov- Taylor solution (jGranot and Sari . 
20021: ISari and Piranl.ll99,5t) . For an extensive review on the GRB blast wave 

ph ysics, seelPiraiiTf 199£ ]^ 

Rees and MeszaroT i 19921 ) and Meszaros and Reei ( 19931 ) suggested syn- 
chrotron radiation as the main radiation mechanism. The synchrotron radi- 
ation assumption naturally implied that the spectrum would soften and fade 
as a power law in ti me, and that an optical and radio afterglo w would be 
present at later times ( Paczvhski and Rhoadsl 19931: Kat d . 1994). The radia- 
tion from both internal and external shocks is fairly well fitted by synchrotron 
and inverse Compton radiation from a high-energy non-thermal electron pop- 
ulation in a strong magnetic field. Below (section 11.31 and II. 4j) I discuss the 
non-thermal acceleration and the possible origin of the magnetic field. 




e [day3 since GRB 99051D] 



Figure 1.2: The jet break of GRB 990510. The break is interpreted as the hmit 
where the relativistic jet is decelerated enough so that the relativistic beaming 
angle ( l/F) becomes larger than the jet opening angle 9j. When this happens, the 
beaming angle covers an increa singly larger ar e a out side the jet and the temporal 
decay will appear faster. From Harrison et al. ( 19991 ). 



It is now well established that the relativistic gamma-ray burst ejecta are 
coUimated. A coUimated outflow is indeed a more compelling scenario since 
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the required energy release from the progenitor is greatly reduced. If the 
emission were isotropic, the cosmological distances would imply that some 
bursts emit more than one solar mass in gamma rays. Such energies are 
hard to produce instantaneously in any stellar model. An observational fact 
that supports the coUimation model is an achromatic break in the power 
law slope of the light curves. For inany afterglows this happen s days to 
weeks after the burst ( Kulkarni et al. . 19991: Harrison et al. . 1999f ) (see Fig. 
II. 2|) . Such a break was suggested and interpreted as the limit where the 
relativistic jet is decelerated enough so that the relati yistic beaming angle 
( 1/T) becomes larger than the iet opening ang le 9j ( Rhoad j . h997[[T99i 



Panaitescu and Meszarosl . Il999l : Sari et all . Il999[ ). When this happens, the 
beaming angle covers an increasingly larger area outside the jet and the 
temporal decay will appear faster. 

One of the big unanswered questions concerns the jet structure. Cal- 
culating the energy budget for a burst requires crucial knowledge of the 
angular shape of the jet. The simplest structure is a jet where the internal 
energy, density and bulk Lorentz factor are constant throughout the jet cone 
( RhoadsL Il999[ ISari et all Il999l ). A more advanced model is the universal 
structured jet where the jet param eters vary smoothly with the angl e mea- 



sured from the iet symmet ry axis (jLipunov et al.L 12001 



Zhang and Meszaros 



rmmetr 
200i). 



Simulations by IZhang et al. 



Rossi et all . 120021 : 



(|2004l ) of the jet 



break-out from a massive Wolf-Rayet star show that the jet consists of at 
least two components (a highly relativistic thin jet and a less relativistic 
cocoon jet). 
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Figure 1.3: Correcting the to tal emitt ed gamma-ray energ y corrected for beaming 
angle. From lFrail et al.l (|2nnih left and lBloom et al.l (j2nn.^ ) (right). 
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Correcting for the jet geometry, the total required burst energy drops 
from 10^^ J to around 10^^ J, not too far from the supernova output. More 
remarkably, in both the uniform structured jet and the universal structured 
jet model, the gamma-r ay energy re l eases for many bursts are narrowly clus- 
tered around 5 x 10^^ j dprail et allboniUBloom et alibnoat) (see FiEf.Oll. 
In the uniform structured jet model, the different light curve break times are 
explained by different opening angles. The universal structured jet explains 
the break time spread by differences in the angle between the line of sig ht 
and the jet symmetry axis ( Rossi et al. . 20021 : Zhang and Meszarosl . 2002[ l. 




1Q48 1 049 1 050 1 051 1 052 1 055 1 Q^* 10^^ 

Energy [erg] 



Figure 1.4: The lAmati et al.l (|2002l ^ and lOhirlanda et al.l (|2004h relations. The 
Amati relation links the total isotropic equivalent gamma-ray energy to the peak 
energy. The Ghirlanda relation {solid line) links the beaming corrected total en- 
ergy to the vF,, peak energy in the cosmological rest frame of the burst. From 



Ghirlanda et al 



{ )eak ei 
200i)- 



Another striking correlation is the Amati correlation. Amati et al. ( 2002t ) 



found a correlation between the peak in uFi, and the total isotropic equivalent 
burst energy. An even tighter correlation was found between the typical pho- 
ton energy and the bea ming corrected gamma-ray output during the burst 
( Ghirlanda et al. . 20041 ) (see Fig. II. 4|) . These correlations have so far been 
purely empi rical, with no viable physical explanation. Recently however, 
Rvdel ( 2005 ) suggested a hybrid model consisting of a strong thermal com- 



ponen t accompanied by a non-thermal component of similar strength. Rvdel 
(200i) made fits of 25 strong bursts and compared the hybrid model to the 
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commonly used Band function ( Band et al. . 1993[ ). The result showed almost 
equally good fits for the two models except for ten of the burst where the 
hybrid model is marginally better. 

In Ryde's hybrid model, the peak of the burst is determined primarily by 
the temperature and is less sensitive on F. In this case, the Amati/Ghirlanda 
correlations have a natural explanation, sinc e for a thermal emit t er the lumi- 
nosity and the temperature are correlated. Rees and Meszaroi ()2004l ) show 
that a correlation close to the observed one arises naturally under certain 
assumptions. 

Thus, it would seem that GRBs are becoming excellent standard candles 
for probing the dark ages of the universe. 



1.3 Particle acceleration in collisionless shocks 



One of the key ingredients in generating what is believed to be non-thermal 
synchrotron radiation from GRB afterglows is a non-thermal, high-energetic 
electron population. Placing an ensemble of electrons with a power-law en- 
ergy distribution function dN{E) oc E'^dE in a homogenous magnetic field 
will result in a synch rotron radiation spectrum wi th a power-law segmeiit 
F{u) oc z/-(P-^)/2 ^g_g_ iLandau and Lifshit jll97,^ and lRvbicki and Li ghtmanlll979l) . 

Many afterglow and Cosmic Ray models assume that elect rons are acce l- 
erated in collisionless shocks by diffusive Fermi acceleration (Fermi, 1949( ). 
In first-order Fermi acceleration, the particles are accelerated as they repeat- 
edly cross the shock transition jump. In the shock rest frame, the incoming 
upstream particles are stochastically deflected by a magnetic field as they 
pass the shock region. In this process, a fraction of the particles are kicked 
to higher energy. These high-energy particles then run into the upstream 
region and a fraction of the particles are reflected into the shock again for 
further acceleration. This iterative process continues in a competitive game 
between energy gain and escape of particles. 

The theory of diffusive shock acceleration predicts that for non-relativistic 
shocks, the resulting particle distribution function converges to a power-law 
with the slope 

Vu + 2vd 3vu ,^ 
p = or s = (LI) 



Vd 



Vd 



fe.g. lAxford et ahlll 9771 iBdllll 9781 and iBIandford and Ostriker] 119781 ). Here 
Vu and Vd are the upstream and downstream plasma bulk velocities. The 
p + 2 is more common in the literature of particle acceleration 

oc 



notation s 



{d^N{p)/dp^ (X p where p is the particle momenta). 
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In the relativistic limit, the Fermi acceleration formalism becomes more 
com plicated. The non-re l ativis tic derivation is based on the diffusion equa- 
tion (jKirk and Schneideii . ll987l ). derived under the assumption that the par- 
ticle distribution function is approximately isotropic in the local plasma 
rest frame. But for relativistic shocks, the bulk flow is comparable to the 
particle velocities and then the particle angular distribution function be- 
comes highly anisotropic near the shock. In this c ase, strong magne t ic fluc - 
tuations downstream of the shock are 

Kirk and Schneider (1987) and Heavens and Drurvl ()l988l ) investigated the 
relativistic problem and found the distribution slopes 



P 



R + 2 
R-1 



3R 



or 



R 



:i.2) 



where i? is the sh ock compression factor. Both analytical dKirk et al. . 200ol ) 
and Monte-Carlo ( Bednarz and Ostrowskil . 19981 : Achterberg et al.r200l[ ) find 
that p converges at 2.23 for F oo (F is here the bulk Lorentz factor). 

The major force of relativistic first-order Fermi acceleration is that it 
predicts indices very close to the ones inferred from observ ations. Estimates 
from a number of GRB afterglows yield p = 2.2 ± 0.2 ( Waxman . 1997bl: 



Berger et al. . 2003[ ). Good agreement and predictions are, however, not the 



same as a scientific proof. Afterg lows exist that have a much larger variety 



in p. E.g. ICampana et al.l (120051 1 find p = 1.3 in the very early afterglow. 



Here I emphasize some of the problems that the Fermi acceleration scenario 
in GRB afterglows is still facing: 

Problem 1 It is very important to stress that Fermi acceleration in collision- 
less shocks is still not understood from first principles. The foundation 
of the acceleration mechanism is based on the test-particle approxi- 
mation. It is assumed that the particles scatter on electromagnetic 
waves but the model does not self-consistently account for the gener- 
ation of these waves. Nor does it account for the back-reaction that 
the high-energy particle distribution have on the electromagnetic field. 
Acceleration that results from currents and charge separation near the 
shock must be probed with a full kinetic approach (e.g. particle-in-cell 
codes). See section IT. 3. II and Chapter |S1 

Problem 2 In all derivations of the relativistic Fermi acceleration mecha- 
nism, the downstream magnetic field is required to b e stronglv turbu- 

lent on scales smaller than the typical gyro-radius (e.g. lOstrowski and Bednar7i2002h. 
This is, however, in conflict with the GRB afterglow synchrotron in- 
terpretation where the high-energy particles are expected to gyrate 
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in circular orbits in a magnetic field with variation scale length much 
longer that the gyro-radii. Either the ansatz of strong downstream tur- 
bulence must be relaxed, with the result that the high-energy particle 
distributi on function is not nearly as universal and possible not power- 
law at all (lOstrowski and Bednar3 . 20021: Niemiec and Ostrowski 2004 ; 
Baringil200,C Or, the radia tion model must be altered to include jitter 



radiation (Medvedev, 2000| ). See also Chapter [7| 



Problem 3 Relativistic Fermi acceleration requires a pre-acceleration mech- 
anism that injects electrons into the iterative acceleration process. The 
pre-acceleration mechanism is not well-known. How large a fraction of 
the electrons that are accelerated greatly affect s our estimates of the 
total GRB energy (jEichler and Waxman . 2005). We defiri e this frac- 
tion as /. Moreover, according to lBaring and Brabvi (j2004^ . agreeable 
synchrotron and/or inverse Compton fits are only attainable when the 
electron population has a significant non-thermal component. This is in 
disagreement with the Fermi process where electrons are injected from 
a dominant thermal pool. The lack of a dominant thermal pool also 
raises a question of how the electromagnetic turbulence is sustained in 
the shock-region. 



Problem 4 This problem is connected to problem 3. In the closest and best 
studied mildly relativistic shock in the Crab Nebula, most of the elec- 
trons r adiate below the expected in jection energy and this means that 
/ ^ 1 ( Eichler and Waxman . 2005| ). The "low" e nergy electrons have 
a pow er-law distribution spectrum 1.3 > p > 1.1 (jWeiler and Panagia . 
1978[ ). This is much lower than what is expected from test particle 
simulations. The high-energy electrons, however, are more consistent 
with slope expected from first order Fermi acceleration p ^ 2.2. 



Problem 5 If the standard Fermi diffusive shock acceleration theory is cor- 
rect, one expects an X-ray halo around the shock. This is because 
higher energy electrons are expected to diffuse further ahead of the 
shock, so the halo should become more extensive at X-ray wavelengths. 
Long eTZI (|2003l ) have investigated high resolution Chandra images of 
the close by Supernova remnant SN 1006. They fail to detect such a 
halo. Instead they see a sharp jump in emissivity at the shock. They 
conclude that either Fermi acceleration is absent in this shock, or some 
kind of shock instability must be operating in shocks that can create 
or amplify a magnetic field with a factor significantly larger than that 
given by the fluid compression, resulting in greater contrast between 
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upstream and downstream emission (|Long et al.l . l2003l ). 



1.3.1 Particle acceleration in PIC simulations 

Clearly, further advances in our understanding of particle acceleration in col- 
lisionless shocks require a full, self-consistent kinetic treatment. This may 
be provided by particle- in-cell (PIC) simulations (see Chapter |2)). Unfortu- 
nately, PIC code simulations are very computationally demanding. There- 
fore, many PIC simulations up to date are one-dimensional. Here, I briefly 
review the state of particle acceleration in PIC code simulations. 

One of the first reports of non-thermal accelerati on in PIC simulation s 



of relativistic coUisionless plasma shocks came from Hoshino et al. ( 1992| ). 



In their (one-dimensional) simulations, plasma, carrying a magnetic field, is 
injected at the leftmost boundary. At the rightmost boundary, the plasma 
flow is reflected and thus collides with itself. In a plasma consisting purely 
of electron/positron pairs, they find that the acceleration is mainly thermal. 
When protons are present the positrons are strongly accelerated. The accel- 
eration is driven by resonant absorbtion o f magnetosonic waves , excited by 
energy dissipated from the gyrating ions. iHoshino et all ()l992h found that 
the positrons could be accelerated to a power-law distribution with slope 
p = 1.5 — 2.5 and that the spectrum extended up to TrriiC^. The reason why 
only positrons and not electrons are participating in the acceleration has to 
do with the polarization of the magnetosonic waves. 

Another mechanism was examined wit h PIC simulations by 
2OOOLI2OO4II . Based on theoretical work bv lOaleev a7a1 (I199,^ , 



Dieckmann et al 



Dieckmann et al. 



2OOOI |200J) used PIC simulations of counter-streaming proton beams on a 



cold plasma background in a transverse magnetic field. They found traces 
of non-thermal particle acceleration up to mildly relativistic energies, offer- 
ing an explanation to the injection problem. One should note, however, that 
their simulations were limited to one spatial dimension and the setup in itself 
is based on many assumptions about the initial conditions. It is not clear 
if the ion beam is injected on a quasi-neutral plasma background. If this is 
indeed the case, there is an excess of correlated positive charges in the sim- 
ulations. This might trigger unrealistic instabilities. The size and duration 
of the simulations are rather limited and it would be interesting to see the 
same simulations carried out in three- dimensions. 

A very well studied a cceleration mechanis r n is t he so called surfatron 
( Katsouleas and Dawson . 19831 : Dawson et al. . 1983[ ). The surfatron is a 
mechanism in which a particle is accelerated while it is trapped by a prop- 
agating large-scale perpendicular electrostatic wave. The waves are driven 
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via the Buneman inst ability tha t arises whe n electron and ion beams drift 
at different velocities ( Buneman . 19581 1959f ). The acceleration of particles 
by trapping in electrostat i c wav es ha s been simulated with PIC sim ulations 
by Hoshino and Shimada ( 2002[ ) and Shimada and Hoshino (12004 ) and has 
been linked to phase space holes or loopholes (jSchmitz et al. 

EM- 



2002) (see Fig. 




Figure 1.5: Phase space plot of ions {top panel) and electrons {bottom panel). 
The sulfation mechanism accelerate electrons by trapping in nonlinear wave modes 
linked to the existence of solit ary electron phase sp ace holes. The result is based 

(|2nn2l ). 



on ID PIC simulations. From Schmitz et al 



It should be noted that several independent simulations have indicated 
that the shock-fronts of collisionless plasma shocks, propagating in an am- 
bient magnetic field, show great time variability (e.g. iLembege et al.l (j2004^ 
and references therein). These simulations include both PIC simulations and 
hybrid simulations. In the latter, the ions are treated as particles an d the elec- 
trons as a massless fluid. In one-dimensional PIC simulations bv iLee et al 



( 2004^ with parameters aimed at supernova shocks, the shock structure was 
foun d to be cyc l ically reforming on ion cyclotron timescales. In the reforma- 



tion, iLee et al.l (|200J) found both electron loophole acceleration but also a 
suprathermal population of ions that could potentially explain the injection 
of high-energy ions into a Fermi acceleration scenario. The acceleration of 
electrons is interesting in the GRB context, and the acceleration of ions is 
interesting in regard to ultra high-energy cosmic rays. 

All the simulations reviewed above have two things in common, 1) they 
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are one-dimensional and 2) they assume that a rather strong transverse mag- 
netic field is present. It would indeed be very interesting to explore whether 
effects such as electron trapping by electrostatic solitary waves will survive in 
three-dimensions, or if other effects become dominant, rendering the previous 
results artifacts of the one- dimensionality. 



Specifically, I would like to express my reservation concerning the surfa- 
tron and loophole acceleration in relativistic collisionless shocks for the fol- 
lowing reasons. The main driving mechanism is the excitation of electrostatic 
waves as a result of the Buneman two-stream instability. However, at rela- 
tivistic velocities, the Weibel two-stream instability has a larger gr owth rate 
than the Buneman instability and will dominate the shock region ( Califanol . 
But the Weibel two-stream instability cann ot be represented cor- 



rectly in simulations with only one spatial dimension. Shimada and Hoshinol 
( 2004^ have shown that the generation of loopholes is suppressed when the 
plasma frequency to cyclotron frequency ratio {upe/uce oc y/n/B) is less than 
10 {upe = (ng^/mgeo)^^^ is the electro n plasma frequency and 0;^^ = qB/nip, 
is the electron cyclotron frequency). iHededal and Nishikawal (|20o3) found 
in three-dimensional simulations that for ujpe/ujce > 10, the initial ordered 
ambient magnetic field becomes curled and even loca lly reversed because of 
the Weibel instability. Hededal and Nishikawal ()2nn5h did find non-thermal 
acceleration, but for other reasons (see Chapter EJ. So for ujpe/uice < 10 the 
solitary electrostatic surfatron acceleration mechanism is suppressed and for 
^pe/^ce > 10 the Weibel two-stream instabili ty is present, which will distort 
the electrostatic wave generation. Finallv. iHede'dal and ^^^I^ B 
found that in the interstellar medium where Upe/uce — 1000, the Weibel 
two-stream insta bilitv evolves unhindered , with no signs of loopholes or sur- 
fatrons (although lHededal and Nishikawa ( 20051 ) do point out the importance 
of larger simulations) . 



Clearly, the time has come for 2D or even 3D simulations to provide a 
more self-consistent explanation for the origin of the ambient magnetic field 
and non-thermal particle acceleration. Evidence from 3D PIC-simulations is 

gathering that suggests that particle acceleration and magnetic fiel d genera- 

tion are two highly connected features of coUisiqnless plasma shocks (IPrederiksen et al 



2002; 



2005 



Silva et al.l. [20031: IPrederiksen et alLbooilNishikawa et al.Ll2003[ iNishikawa et al.l 



Hededal et all booi iHededal and Nishikawal . 1200,^ . I save the discus- 



sion of these 3D PIC simulations of particle acceleration and the connection 
to field generation to section 11.4.11 below. 
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1.4 Magnetic fields in gamma-ray bursts 



A second crucial ingredient in generation of radiation in GRB afterglows is 
the presence of a st rong magnetic field. For a review on the role of magnetic 



fields in GRBs, see iPiranI f)20n5aD . 



The general working assumption is that the energy that resides in the 
magnetic field and in the non-thermal electrons may be parameterized by 
the equipartition parameters and eg, and that the non-thermal electrons 
follow a power-law distribution with slope p (where and eg are defined 
as the fractions of the total internal energy of the shock that are deposited 
in magnetic energy and kinetic energy of the electrons, respectively). It is 
generally assumed that these parameters are constant through the shock and 
even throughout the duration of the afterglow. The values are observationally 
determined by localizing certain characteristic break frequencies in the spec- 
tra. From low to high frequencies these are the synchrotron self-absorbtion 
frequency, the sync hrotron frequency of the typica l electro n, and the self ab- 
sorption frequency ( Sari et al. . 19981 : Piran . 19991 2005a 3)- Regarding the 
mangetic field, the typical val ue of the equipartition parameter in the af- 
terglow is 6r = 0.0001 - 0.1 (IWaxmanI Il997bl : IWiiers and Galamal . Il999l : 
Panaitescu and Kumar . 20021 : lYost et al.L [2OO3I ). This value may be trans- 
lated to a magnetic field of the order of lO"'^ T (1 G) in the afterglow shock. 
In the interstellar medium the typical magnetic field strength is of the order a 
few 10~^° T (few /iG). According to the relativistic Rankine-Hugoniot plasma 
shock jump conditions (lTaubl .1 Il948f ). shock compression can only give a fac- 
tor of 4:r shock and t his is clearly far from enough to match the interpretations 



1999h . This leaves us with two 



from observations (jGruzinov and Waxman . 
possibilities for the origin of the magnetic field. One possibility is that the 
magnetic field is generated or amplified in the shock by microphysical insta- 
bilities and one is that the magnetic field is carried with the outflow from the 
progenitor. A magnetic field that originates from the progenitor and is frozen 
into the ejecta might account for the magnetic field in the internal shocks. 
But as the plasma shell expands, the magn etic field is diluted arid diss ipated 
to well below the auticpated values (e.^. l M.Hv.d.v „.H r,n.lJ I Il999>) . Ad- 
ditionally, there exists a question of how to transport a magnetic from the 
ejecta and into the shocked ISM (all though theoretically it may happen via 
the Rayleigh- Taylor instability). Hence, the magnetic field responsible for 
the afterglow is most likely to be generated in situ in the shock. 

Under collision of two relativistic plasma populations, the particle phase 
space is extremely anisotropic. Naively one could argue that in the ab- 
sence of particle collisions, the two plasma populations would stream right 
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through each other. The anisotropy is, however, unstable to several plasma 
instabilities, including the electrostatic Buneman two-stream instability and 
the electromagnetic Weibel two-stream instability. At rel ativistic shocks, 

the latter has the largest gro A vth rate and will dominate (ICalifanol. 120021: 

Hededal and Nishikawa . 2005| ). Gruzinov and Waxman ( 199£ ) and lMedvedev and Loebl 
(jl999| ) suggested that the Weibel two-stream instability could generate a 
strong magnetic field in the shock region. Since many of the results that 
I present in this thesis are connected to the Weibel two-stream instabil- 
ity, it is worthwhile to explain the nature of t he instability i n some detail. 
The following description i s base d on papers by Weibell ( 1959f). Freid ( 1959t ) , 
Medvedev and Loeb l ( 1999t ). and lWiersma and Achterberg I ( 2004i ). 



Weibell (|l959t ) suggested that if an isotropic plasma population is anisotrop- 



ically perturbed, relaxation will lead to a growing transverse magnetic field, 
even i n the absence of an external electromagnetic field. The same year, iFreid 
( I959I ) gave a physical interpretation, where the anisotropic perturbation was 
described as an actual two-stream configuration. Since there will always be 
infinitesimal magnetic perturbations in a plasma. Fried suggested that defiec- 
tion of the electrons (by the Lorentz force) in such a fiuctuating magnetic field 
will create currents that can amplify the initial magne tic perturbation. Fig- 
ure 11.61 shows a schematic drawing of this mechanism ( Medvedev and Loebl . 
19991 ). In the center of mass rest frame, two oppositely directed electron 
beams collide at x = (with bulk velocities Vi^ = —v^ and V2x = so 
that the net current is zero). The ions are treated as a homogenous non- 
interacting background, to ensure charge neutrality. At x=0, a magnetic 
perturbation is initially present with Bz = Bq cos ky. This perturbation de- 
fiects left streaming and right streaming electrons into anti-parallel currents 
with a distance comparable to the wavelength of the magnetic perturbation. 
According to Ampere's law these currents represent a curl in the magnetic 
field. This curl amplifies the initial perturbation, which, in turn, collects 
even more electrons into the current channels. This positive feedback results 

Ve/C 



L0< 



pe 



i n an instab ility where magnetic fiuctuations grow with a rate a 
(Freid, 19591 ). Again, Upe is the electron plasma frequency defined above. 

The relati vistic generalization is not trivial. The problem have been in- 
vestigated bv lYoon and DavidsonI (jl987l ) who find the maximum growth rate 
to be (J = Wpe/r^/^ where F is the bulk relativistic Lorentz factor. The result 
can be understood intuitively by evaluating the non-relativistic expression in 
the limit where Ve —>■ c and Upe — > {nq'^ /Trrigeo 



1 1/2 



Medvedev and Loebl ( 19991 ) estimated that the instability would saturate 



at ~ 10 ^ — 10 ^ if only the electrons participa te in the instability, and 
< 10^^ if also the ions take part in the instability. IWiersma and Achterberg 
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Figure 1.6: Schematic of the Weibel two-stream instabihty. Counter-streaming 
electrons are collected into opposite directed current channels by a magnetic 
perturbation. The curre nt channels amplify the initial perturbation. From 
Medvedev and Loebl (Il99fll 'l. 



( 2004 ) found from a linear analysis that the two-stream instability in an 
electron-proton plasma shock has an early end where < me/mj and that 
the wavelength of the most efficient mode for magnetic field generation equals 
the electron skin depth. Numerical simulations of the instability have shown 
that this limitation not true for the non-linear stage. 



1.4.1 Magnetic field generation in PIC simulations 

To investigate the non-linear stage of the Weibel two-stream instability, self- 
consistent kinetic particle-in-cell (PIC) simulations are necessary. In Chapter 
121 1 briefly describe the theory behind PIC codes with their advantages and 
disadvantages. 

The ffist PIC simulations with direct foc us on the Weibe l insta bility in 
an astrophysical context were performed by iKazimura et al. I (Il998t) . They 



were interested in the plasma wind interaction in millisecond binary pulsars. 
In two-dimensional simulations of size 6.6 x 106.6 electron skin depths, they 
investigated the collision of two mildly relativistic pair-plasma winds {vq = 
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Figure 1.7: The structure of current channels generated by the Weibel two- 
stream instabiUty under the cohision of two pah'-plasma shells {left panel). The 
corresponding amount of generated inagnet ic end electric field is close to equipar- 
tition (right panel). From Silva et al. ( 2003h . 



The collisio n of pair plasma sh ells was investigated with three-dimensional 
simulations by Silva et al. ( 2003[ ). In a simulation box with size 25.6 x 25.6 x 
10.0 electron skin depths, they explored the collision of pair-plasma shells 
with relativistic Lorentz factors ranging from F ~ 1 to F = 10. They found 
that the generated magnetic field reached maximum of ~ 0.1 and after- 
ward relaxed to ~ 0.0001 — 0.01 (see Fig. II. 7|) . I emphasize that in these 
simulations the boundaries in the flow direction are periodic. It is not clear 
what this implies, but it can potentially introduce non-physical, stabilizing 
fee d-back, since each curr ent-fi lament is feeding itself. 



Nishikawa et al.l ()200,'^ and lNishikawa et al.l ()2005r i have performed three- 



dimensional simulations of both pair-plasma and electron-proton shocks in a 
numerical box with size 17.7 x 17.7 x 66.7 electron skin depths (corresponding 
to 4 X 4 X 14.9 ion skin depths for ion-electron mass ratio mj/mg = 20) (see 
Fig. ll.8p . With these simulations they were able to follow primarily the linear 
stage of the instability at the front of the shock ramp (the ambient and jet 
electron populations are not the rmalized to a sin gle population in the sim- 



ulations) . W ith the s imulat i ons iNishikawa et a .1 confirmed the growth rat e 



predicted by IWeibell ()l 959l ). iFreidI ()l959l ) and iMedvedev a,nd Loebl ()l 999l ). 

They also found signs of electron acceleration connected to the instability, 
but whether the nature of this acceleration is non-thermal or merely a ther- 
malization is not clear. 
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Figure 1.8: The electron density {left panel) and current density (right panel). 
The arrows show the induced magnetic field. From Nishikawa et al. (|2005l ) . 
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Figure 1.9: The left hand side panel shows the longitudinal electron current 
density through a transverse cut early in the shock. The right hand side panel 
shows the ion curren t deeper in the shock. T he arrows represent the transverse 
magnetic field. From Frederiksen et al. (12? ' 



Frederiksen et a1.l \2m± 12004 ) and iHededal et all ()2004 ) have investi 



gated the non-linear evolution the two-stream instability in electron-proton 
shocks. Details are givei i in Chapter El an d Chap ter IHl of this thesis. Here 
I give a short summary. iFrederiksen et al.1 (j200J) performed simulations of 
electron-proton shocks with F = 3. The size of the simulation box was 
40 X 40 X 160 electron skin depths (corresponding to 10 x 10 x 40 ion skin 
depths for ion-electron mass ratio mj/mg = 16). The results of these sim- 
ulations showed how Weibel generated ion current filaments were collected 
into increasingly larger patterns in the non-linear stage. This phenomenon 
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has been investig a ted both analytica l ly and numerically in great details by 



MedvedevI f)2005[ ). iFrederiksen et al.1 (j200J) also found that the ion current 
filaments are partially Debye shielded by shock-heated electrons. The elec- 
trons thus act as insulators for the ion-current filaments, making these rather 
robust (see Fig. II. 9|) . The results showed that a magnetic field was generated 
{sb ~ 0.0 5) and sus t ained for many ion skin depths. This answers a concern 
raised bv iGruzino 3 (Il999> ). who estimated that the magnetic field can only 
be sustained for roughly one electron skin depth. 
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Figure 1.10: Tlie high-energy electrons (red colored dots) are found whe re the 
ion-current filaments are strongest (blue colors). From Hededal et al. ( 20041 ). 



Hededal et al.l 1)2004^ found that not only can the Weibel two-stream in- 



stability account the generation of a strong magnetic field, but it appears 
that non-thermal electron acceleration is a natural consequence of the pro- 
cess. They found that high-energy electrons are spatially connected to the 
ion current filaments (see Fig. ll.lOl and Chapter EI). The acceleration and 
deceleration of electrons is local and instantaneous. This is in contrast to the 
recursive process of Fermi acceleration. 



1.5 Summary and Thesis Outline 

The radiation from ORB afterglows is produced in relativistic, collisionless 
plasma shocks by two key ingredients, namely 1) a population of highly en- 
ergetic, non-thermal electrons and 2) a strong electromagnetic field. All our 
knowledge of GRBs are based on this radiation. Nevertheless, the questions 
of how the electrons are accelerated, and what the exact generation mech- 
anism and nature of the electromagnetic field in the shock is, have not yet 
been answered in a self-consistent way and are still open questions. 
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Even though magnetic field generation by the Weibel two-stream insta- 
bility seems necessary (and also unavoidable) in collisionless shocks in GRB 
afterglows, it have not yet been possible to verify or discard it from observa- 
tions. One way to investigate on the magnetic morphology in these shocks 
is to measure the pol a.rization of aft e rglow s, although this also relies heavily 



on the jet structure. iLazzati et al.1 (j2004^ found that polarization of GRB 



020813 is not very well fitted with a homogenous jet with shock generated 
field. I stress that to draw any conclusions from observations about the mag- 
netic field and particle acceleration, it is vital to have a firm understanding 
of collisionless shocks, and the generation of radiation in these. 

It is the goal of this thesis to shed light on the microphysics of collisionless 
plasma shocks, mainly in the context of gamma-ray burst afterglows. Using 
a particle-in-cell code that works from first principles, the aim is to obtain 
insight in, and explain the origin of, the magnetic field in these shocks, the 
nature of the field, and how it may infiuence the radiation emitted from GRB 
afterglows. 

The thesis is divided into M chapters: 



• In Chapter |21 1 describe the particle-in-cell code that have been used 
and how radiative cooling has been implemented. 

• Chapter|21presents the results of simulations of the non-linear evolution 
of the Weibel two-stream instability. This chapter is based on the paper 
by Frederiksen, Hededal, Haugb0lle and Nordlund (2004). 

• In Chapter ^ I compare two- and three-dimensional simulations, and 
present results of large-scale shock simulations. This work has not yet 
been submitted to a scientific journal. 

• Chapter El investigates the microphysics of collisionless plasma shocks 
in the presen ce of an ambient magnet i c field . This chapter is based on 
the paper by iHededal and Nishikawa (1200,^ . 

• In Chapter IHl I present a new particle acceleration mechanism that 
differs from Fermi acceleration. This paper is based on the paper by 
Hededal, Haugb0lle, Frederiksen and Nordlund (2004). 

• In Chapter [7| I describe the development and test of a new and powerful 
numerical tool, which may be used to create radiation spectra directly 
from PIC simulations. This work has not yet been submitted to a 
scientific journal. 
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• Chapter |S1 describes the foundations of a next generation particle-in- 
cell code that includes photons and various scattering processes. The 
chapter includes some preliminary test results with Compton scatter- 
ing. 

• In Chapter 121 1 collect all the pieces from the thesis in a summary and 
conclusions. Here I also discuss the future of the line of work that I 
have presented in this thesis. 



Chapter 2 

The Particle— In— Cell code 



In this Chapter I briefly describe the kinetic particle-in-cell code, and justify 
why it is important compared to a fluid description. I also discuss some 
of the hmitations of a kinetic numerical description. Finally, I describe in 
some detail the derivation, implementation and testing of a radiative cooling 
mechanism in the code. 

2.1 Kinetic or fluid description of collision- 
less shocks. 

The mean free path for a 90° Coulomb deflection of an electron moving with 
relativistic momentum Tvefrie in a plasma (with density n) is 



See Appendix^ for details. For a relativistic GRB jet that expands into the 
interstellar medium (ISM), we may use this expression to estimate the typical 
mean free path for Coulomb collisions between the jet and ISM particles. 
With an ISM density n ^ 10^m~^ and a jet bulk Lorentz factor 7 = 5 (f ~ c), 
the mean free path for Coulomb collisions is lO^^m. This is of the order of a 
billion times the expected size of the flreball. Thus one might naively expect 
a relativistic jet to expand unhindered through the ISM. This is however, 
in direct contrast with observations, where gamma-ray burst afterglows are 
described by synchrotron emission from decelerating relativistic shells that 
collides with, and heats, an external medium. Lack of interactions would 
pose serious problems in explaining the particle acceleration and origin of 
the magnetic fleld, which is needed to produce the observed synchrotron 
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radiation (e.g. Waxman 1997a and Sari et al. 1998| ). Not surprisingly, the 
coUisional interaction agent must be found in the microphysical p rocesses at 
play between particles and electromagnetic fields ( Sagdeevl . Il966t ). 

From this discussion it is clear that a treatment of the jet/ISM interaction 
needs to be established. For this purpose we need a theoretical framework. 
The use of the Magneto-Hydro-Dynamic equations (MHD) in this context is 
discarded by several arguments: 

• The low collision rate cannot provide the equilibration of ion and elec- 
tron energies sufficiently fast for the plasma to behave as a fiuid. This 
is the case even f or the "low energy" shocks associa ted with super- 



19931 : IVinkl . 120041 ). Observations 



nova remnants ([Draine and McKe^ 
are consistent with an energetically important population of acceler- 
ated part i cles s uperimposed on a low energy background population 
Gruzinovl (I2OOIL 



MHD shocks are stable and do not generate magnetic fields. In MHD 
shocks, magnetic fields are only compressed, with a resulting field 
strength that is orders of magnitudes smalle r than what is required 
by the synchrotron model of GRB afterglows (jGruzinovl . 120011 ). 



2.2 From first principles 

The solutions that we are seeking fall into the kinetic and highly non-linear 
regime of plasma physics. In this case, we need to work from first principles, 
by solving the Maxwell equations with source terms for the electromagnetic 
fields, together with the relativistic equation of motion for charged particles 



^ ^ dE 
V X B 
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and 

=g(E + v X B). (2.3) 

at 

Here eg and /iq are the constants of electric permittivity and magnetic per- 
meability of vacuum, with c^/ioCo = 1- m and q are the mass and charge of 
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a particle of a given species, v is the velocity vector and 7 = (1 — t>^/c^)~^/^ 
is the relativistic Lorentz factor. The source terms, J and p, in the Maxwell 
equations are determined by the particles in the simulations. 

We wish to find a general solution to the coupled differential equations of 
Eg. 12.21 and Eg. 12.31 For a given set of initial/boundary conditions with, say, 
10^^ particles, this is not analytically possible. Numerically, however, solving 
a scaled-down ve rsion of the same problem is doable, with particle-in-cell 
(PIC) codes (e.g. Birdsall and Langdon ^^). Much like in a real plasma, a 



PIC code integrates the trajectories of a large number of charged particles in 
both external and self-induced electromagnetic fields. Some limitations exist 
in this approach. Some of the main differences between a PIC simulated 
plasma and a real plasma are: 

• The number densities of real space are many orders of magnitudes 
larger than what can be fitted in a computer: A 10 x 10 x 10 m^ cube 
of the ISM contains approximately 10^ charged particles, and even this 
is barely computational affordable today. Therefore, in the simulations, 
each particle is a macro-particle that represents a large number of real- 
plasma charges. Each macro-particle keeps the same charge to mass 
ratio as the individual particles it is made of. 

• Even though continuous in space and momentum space, the particles 
positions are discretized in time. 

• The electromagnetic fields are discretized is space as well as time. The 
Maxwell equations are integrated on a fixed numerical grid and the 
interactions with the particles in a given grid cell are done via interpo- 
lations from grid to particle positions and vice versa (hence the name 
particle-in-cell). The electromagnetic field components a^d source- 
terms are staggered and distributed on a 3D Yee lattice (|^, Il96fil ). 



This gives a resolution improvement that corresponds to a factor 16 in 
computing time fFig. 12. 1|) . 

Many plasma processes evolve on time scales that are proportional to 
the plasma frequency r oc u"^ and on length scales that are propor- 
tional to the skin depth 6 = c/ujp. Therefore, a large spatial and tempo- 
ral span exists in plasma processes that are dominated by respectively 
ions and electrons. To comply with the limitations in computational 
resources it is convenient to compress the dynamical ranges by reduc- 
ing the ion (proton) to electron mass ratio mj/mg from the real value 
(1836) to 15-30. This is clearly an approximation but we have per- 
formed tests that have shown that the results show good convergence 
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Figure 2.1: The staggered mesh. In a grid-ceh, the particle densities are cell 
centered, the electric field and current source terms are face centered and the 
magnetic field components are edge centered. Distributing the field variables on a 
staggered mesh gives a factor 2 increase in resolution. In 3+1 dimensions this is a 
factor 16 in computing time. 



for mass ratios above 15-30. For lower mass ratio, the results are still 
qualitatively correct. 

• The maximum temporal and spatial scale-lengths in PIC simulations 
are limited because it is important to resolve microphysical plasma os- 
cillations. Since the electron plasma frequency Upe is often the limiting 
factor, we normalize time with respect to the oscillation period uo~^ and 
the space with respect to the electron skin depth c/wpe. The plasma 
frequency is defined as ujpe = (ngg^/meeo)"^^^, and thus the plasma 
density rig determines the re-scaling. 

Despite these constraints, the PIC code representation of a plasma is far 
more fundamental than the MHD approximation. Still, PIC simulations are 
computationally challenging and fully three-dimensional experiments have 
only become practically possible within the last few years. 

The PIC-code implementation I use is based on a non-relativistic code 
developed by Dr. Michael Hesse. The code was initially develo ped for simu- 
lating reconnection topologies in the context o f space weather dHesse et al. . 
1999l li. The code was later made relativistic by Frederiksen ( 2nn2[ l as part of 
a masters project. Since then it has been use d mainly for numerical plasma 



shock experirnents related to GRB afterglows ( Frederiksen et al.L 120021 12004 ; 
Hededal et al.L 120041 iHededal and Nishikawal . 120051 ). As part of the current 
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PhD project, radiative cooling has recently been included to the code. This 
is important for investigation of particle acceleration and generation of radia- 
tive spectra. A new PIC code, which includes photons and several scattering 
mechanisms, is under development (see Chapter |H)). 



2.3 Cooling of an accelerated charge 

In this section I derive and describe the implementation of radiative cooling 
in the PIC-code. Radiative cooling is essential for highly relativistic parti- 
cle dynamics and especially for experiments aimed at investigating particle 
acceleration. Before I describe the derivation and implementation an expres- 
sion for the energy radiated from an accelerated charged particle into the 
PIC-code, I briefly revive the concept of retarded time and space. 



2.3.1 Retarded time and position 

Let a particle be in the position ro(t) at time t fFig. 12. 2j) . At the same time, 
we observe the electric field from the particle in the position r. However, 
because of the finite propagation velocity of light, we observe the particle at 
an earlier position ro(t') where it was at the retarded time t' = t — St' = 
t — R(t')/c. Here -R(t') = |r — ^o(t')\ is the distance from the charge (at the 
retarded time t') to the observer point. 




Figure 2.2: Definition of the retardation of a particles position. From an observers 
point r we see the particle at the position ro(t') where it was at the retarded time. 
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2.3.2 Radiated power 

To skip a trivi al, but rather long derivation, we adopt the expression from 
Jacksonl (|l999l ) for the retarded electric field from a charged particle moving 
with instant velocity (3 under acceleration 0, 
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velocity field 



acceleration field 



Here, n = R(t')/ l^-l^OI is a unit vector that points from the particles re- 
tarded position towards the observer. The first term on the right hand side, 
containing the velocity field, is the Coulomb field from a charge moving with- 
out influence from external forces. The second term is a correction term that 
arises if the charge is subject to acceleration. Since the velocity-dependent 
field is falling of as while the acceleration-dependent field falls off as R~^, 
the latter becomes dominant when observing the charge at large distances 
(i? 3> 1). This term is therefore often referred to as the radiation term. The 
corresponding magnetic field is given by 
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The energy W per unit area dA per unit time dt that is radiated from 
the accelerated particle is given by the Foynting flux S 

|2 



d^W E X B 
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dtdA 
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(2.6) 



from which we define the energy per unit time, received through a unit solid 
angle element dQ about n 
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(2.7) 



A note of caution must be added here. The Poynting vector is related to the 
observer time t but we are interested in the radiated power measured at the 
particle's retarded time t'. Thus, to get the total emitted power we multiply 
with the correction term dt/dt' (see Chapter [3 for details) 
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Rather tedious vector algebra and integration over all directions n through 
the solid angle dQ gives us the total power radiated by the particle in a time 
interval dt' (see Appendix |Bl for details) 



rad 



H — ^ 
ovr 

Gtt 



/3P 



76/32 (l-Z^^sin^^) 



(2.9a) 
(2.9b) 



where 6 is the angle between the particles velocity vector f3 and its accel- 
eration vector f3. When sin 6* = 1 <(=^ /3 _L /3, we recognize the solution 
for synchrotron radiation (magnetic bremsstrahlung) . In the other limit, 
sin 6' = -vv- /3 II /3, we recover the res ult for bremsstrahlung. The results 
above may be found in many textbo oks ( Landau and Lifshitj . ll975[ l iJacksonl . 
I999I : iRvbicki and Lighti^ . Il979> ). 

In Chapter [7| we deal with radiation from relativistic particles in more 
details. 



2.3.3 Implementing the radiative cooling 

Before implementing Eq. 12.91 into the PIC-code, it is fruitful to derive the 
expression for a particles energy as it looses momentum to radiation. 

From energy conservation we demand that the energy radiated from a 
particle must equally correspond to a loss in the particles kinetic energy 
Ekin = mc^{'j - 1) 



77 ^rad 

dt 

d'j{t) /iog^ 



/(/52(l-/32) + |/3-/3r), (2.10) 



dt Gnmc 

which is the relativistic counterpart to the Larmor formula for radiated 
power. In order to continue we need an expression for 0. The only force 
acting on the particle is the Lorentz force (we deal with the radiation reac- 
tion force in Section I2.3.4|) 

^c?(7v) ^ g(E + vxB) (2.11) 
dt 

dv V7^ / dv\ ,^ , , 

^^*+^(^-7ft) = + (2.12) 

Here we have expanded the right hand side and used the relation 
d^ d{l - f3 ■ f3)-'/' 7' dv 

Tt = It = ^^-^^^ 
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Unfortunately, we cannot provide a general solution to Eq. I2.1(JI and Eq. 
12.121 We may, however, simplify the equation system by looking at a particle 
moving with 7(t = 0) = 7o at an angle a to a homogeneous magnetic field 
(E = 0). In this case, the Lorentz force is always perpendicular to the 
velocity vector and the term v-v vanishes. Equation 12. 101 and Eq. 12. 121 then 
reduce to 



di{t) /iog 



4 



dt Gnm'^c 



s\n\a)B\^\t) - 1), (2.14) 



where we have used that v'^{t) = c^(l — 7 ^(t)). This differential equation 
has the formal solution 

7(t) ^ , . (2.15) 



(7o + l)exp Jm^-,sm'{a)BH 



70 



If we in Eq. 12.141 assume that 7 ^ 1 ^ f ~ c, the solution becomes simpler, 
but less valid for low 7 (see fig. 12. 3p 

7(t) = . (2.16) 

In Section l!^.3.4I I describe the implementation of the radiative damping from 
Eq. 12.91 into the PIC-code and find good agreement with Eq. 12.151 

Eq. 12.161 gives us an indication of when momentum losses to radiation 
becomes an important issue. If a physical process of interest occurs over a 
time span T, radiation may be neglected whenever 

kT<1, (2.17) 

where 

-- sm\a)B^^o. (2.18) 



We may also estimate the time ti/2 it takes for a particle to loose half of its 
energy 

mc\j{h/2) - 1) = l/2mc2(7o - 1) 
~ /iog^B27o7o + 2 7o»i /iog^527o' 
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Figure 2.3: The two versions of ^{t) - Eg. 12.151 and Eg. 12.161 - compared. The 
main difference is that in deriving Eg. 12.161 w ~ c was assumed instead of the 
correct v = csjx — in Eg. 12.151 Thus, Eg . 12 . 1 6l b ecomes incorrect as 7(t) — > 1. 
The units on the time axis is given in terms of the initial orbital period. 



2.3.4 Implementing the reaction force 



When a particle emits radiation, it looses both energy and momentum to 
the emitted photons. Thus, in the ory, we need to add an extra force term 
in Eg. 12.121 However, as stated by Jackson ( 19991 ) "a complete satisfactory 
classical treatment of the reactive effects of radiation does not exist" (and it 
is noteworthy that Jackson only discuss the issue in the very last Chapter). 
One approach is the Abraham- Lorentz equation of motion 
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~~dt df^ 



horentzi 



(2.20) 



but this solution is limited to periodic motions and have runaway solutions. 
A consistent, but rather extensive , implementati o n of t he reaction force in a 
PIC code have been described bv lNoeuchi et al.l (j200J). 

Instead we derive a more empirical solution to the problem that is suitable 
for a numerical implementation. We recall the expression for the particles 
kinetic energy lost to radiation 
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Figure 2.4: The trajectory of a relativistic particle in a transverse homogeneous 
magnetic field. The particles motion is integrated by the PIC code. The left panel 
shows the orbit without damping whereas the right panel shows the orbit of the 
particle as it looses momentum to radiation according to Ea. l2.23l Initially 70 = 60, 
K = 0.005 and the orbit is integrated over over the time interval T = [0, 1000]. 



We are looking for a way to modify the particles momentum vector u = 7V 
that we may implement in the PIC-code. We allow ourself to make a discrete 
representation of Eq. 12.211 



7-7 
At 
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rad 
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[2.22] 



where 7 is the particle Lorentz factor from the code and 7 is the radiation 
corrected Lorentz factor. The corresponding radiation corrected momentum 
is |£i| = 0^/7^ + 1. 

We know that in an observers frame, the radiation from a relativis- 
tic particle is concentrate d around the direction of its velocity vector (e.g. 
Landau and Lifshitzl IQTli l. Therefore, we can make the assumption that the 



radiation reaction force is directed opposite to the momentum vector. In this 
case, the change in momentum occurs only along the momentum vector u 
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Figure 2.5: A comparison between cooling times in tlie theoretical approach of 
Eq. l2.15l (rerf lines) and the PIC-code implementation of Eg . [2?23l ( black lines). The 
small discrepancies between the two approaches are mainly caused by numerical 
integration inaccuracy. The time is in arbitrary units. 



and we have 
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|u| |u| 



U = U- — r = U 

u 



(2.23) 



where Prad is given by Eq. 12.91 We have implemented Eq. 12.231 in the PIC 
code. To verify the correctness of the method, we test against the result from 
the analytical approach, Eq. 12.151 We make a simulation run with E = 
and a homogeneous magnetic field B. A single particle is setup with a given 
velocity perpendicular to the magnetic field (sin^ a = 1). The integrated 
path of the particle may be seen in Fig. 12.41 The figure shows that as the 
particle cools to lower gamma it approaches a state where cooling becomes 
negligible and the orbit becomes semi stable [kT <C 1). 

A comparison of the PIC-code cooling rate against Eq. l2.15l mav be found 
Fig. 12.51 Five different combinations of 70 and Bq has been tested. Good 
agreement is found between the two approaches. The minor discrepancies 
that do appear are mainly caused by integrated interpolation errors in the 
second order scheme used in the PIC-code. 



Chapter 3 

The Weibel two-stream 
instabihty 



In this chapter I present results from three-dimensional particle simulations of 
collisionle ss shock formation , with relativistic counter-streaming ion-electron 



plasmas (jFrederiksen et all 12004^ . Particles are followed over many skin 



depths downstream of the shock. Open boundaries allow the experiments 
to be continued for several particle crossing times. The experiments confirm 
the generation of strong magnetic and electric fields by a Weibel-like kinetic 
streaming instability, and demonstrate that the electromagnetic fields prop- 
agate far downstream of the shock. The magnetic fields are predominantly 
transversal, and are associated with merging ion current channels. The to- 
tal magnetic energy grows as the ion channels merge, and as the magnetic 
field patterns propagate down stream. The electron populations are quickly 
thermalized, while the ion populations retain distinct bulk speeds in shielded 
ion channels and thermalize much more slowly. The results help us to re- 
veal processes of importance in collisionless shocks, and may help to explain 
the origin of the magnetic fields responsible for afterglow synchrotron/jitter 
radiation from Gamma-Ray Bursts. 



3.1 Introduction 

The existence of a strong magnetic field in the shocked external medium is 
required in order to explain the obser ved radiation in Gamma- Ray Burst af- 
terglows as synchrotron radiation (e.g. Panaitescu and Kumar . 20021 ). Nearly 



collisionless shocks, with synchrotron- type radiation present, are also com- 
mon in many other astrophysical contexts, such as in super-nova shocks, and 
in jets from active galactic nuclei. At least in the context of Gamma-Ray 
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Burst afterglows the observed synchrotron radiation requires the presence of 
a stronger magnetic field than can easily be explained by just compression 
of a magnetic field already pre sent in the external medium. 

Medvedev and Loebl (Il999> ) showed through a linear kinetic treatment 



how a tw o-stream magnetic instability - a genera lization of the Weibel in- 



Yoon and Davidsonl . |l987|) - can generate a strong 



stability (IWeibell . 1195 

magnetic field (e^, defined as the ratio of magnetic energy to total kinetic 
energy, is 10~^-10 7^ of equipartition val ue) in coUisionless shock fronts 



see 



also discussion in iRossi and R^ . 120031 ). We note in passing that this in- 
stability is w ell-known in other plasma physics discip lines, e.g. laser-plasma 



interactions (jYang et al. . I1992I : ICalifano et all. Il998h. an d has been applied 



in the context of pulsar winds by Kazimura et al.l ( 1998[ ) 

Using three-dimensional particle-in-cell simulations to study relativistic 
coUisionless shocks (where an exte rnal plasma impac t s the s hock region with 
a bulk Lorentz factor F = 5 — lOl. lPrederiksen et al. ( 2002| ). Nishikawa et al. 



( 2003[ ). and Silva et aD (j2003[ ) investigated the generation of magnetic fields 
by the two-stream instability. In these first studies the growth of the trans 
verse scales of the magnetic field was li mited by the small sizes o f the compu- 
tational domains. The durations of the Nishikawa et al. ( 2003[ l e xperiments 



were l ess than particle travel times through the experiments, while ISilva et al. 



ther. iFrederiksen et all (|2002i l and iNishikawa et all (|2003h used electron-i 

(|2003t) 



( 20031) used periodic boundary con ditions in the direction of streaming. Fur- 



-lon 

were 



(e~p) plasmas, while experiments reported upon bv ISilva et al 
done with e~e^ pair plasmas. 

Here, we report on 3D particle-in-cell simulations of relativistically counter- 
streaming e~p plasmas. Open boundaries are used in the streaming direction, 
and experiment durations are several particle crossing times. Our results can 
help to reveal the most important processes in coUisionless shocks, and help 
to explain the observed afterglow synchrotron radiation from Gamma-Ray 
Bursts. We focus on the earliest development in shock formation and field 
generation. Late stages in shock formation will be addressed in successive 
work. 



3.2 Simulations 

Experiments were performed using a self-consistent 3D3V (three spatial and 
three velocity dimensions) electromagnetic parti cle-in-cell code orig inally de- 
veloped for simulating reconnection topologies ( Hesse et al.L Il999[ ) . redevel- 



oped by the present authors to obey special relativity and to be second order 
accurate in both space and time. 
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The code solves Maxwell's equations for the electromagnetic field with 
continuous sour ces, with fi elds and field source terms defined on a staggered 



3D Yee-lattice (jYed . 1966t ). The sources in Maxwell's equations are formed 



by weighted averaging of particle data to the field grid, using quadratic spline 
interpolation. Particle velocities and positions are defined in continuous 
(r, 7v)-space, and particles obey the relativistic equations of motion. 

The grid size used in the main experiment was {x, y, z) = 200 x 200 x 
800, with 25 particles per cell, for a total of 8 x 10^ particles, with ion to 
electron mass ratio rrii/me = 16. To adequately resolve a significant number 
of electron and ion skin-depths (5e and 6i), the box size was chosen such that 
106i ~ 40(5e and Lz ~ 405j ~ 1605e- Varying aspect and mass ratios 



were used in complementary experiments. 




Figure 3.1: The left hand side panel shows the longitudinal electron current 
density through a transverse cut at z = 100, with a small inset showing the ion 
current in the same plane. The right hand side panel shows the ion current at 
z = 600 = 306i, with the small inset now instead showing the electron current. 
The arrows represent the transverse magnetic field. Both panels are from time 



1200 (240u;, 



pe 



Two counter-streaming - initially quasi-neutral and cold - plasma popu- 
lations are simulated. At the two-stream interface (smoothed around z = 80) 
a plasma [z < 80) streaming in the positive z-direction, with a bulk Lorentz 
factor r = 3, hits another plasma {z > 80) at rest in our reference frame. The 
latter plasma is denser than the former by a factor of 3. Experiments have 
been run with both initially sharp and initially smooth transitions, with es- 
sentially the same results. The long simulation time allows the shock to grad- 
ually converge towards self-consistent jump conditions. Periodic boundaries 
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are imposed in the x- and ^/-directions, while the boundaries at 2 = and 
z = 800 are open, with layers absorbing transverse electromagnetic waves. 
Inflow conditions at z = are fixed, with incoming particles supplied at a 
constant rate and with uniform speed. At 2; = 800 there is free outfiow of 
particles. The maximum experiment duration is 480 (where Upe is the 
electron plasma frequency), sufficient for propagating F ^ 3 particles 2.8 
times through the box. 

3.3 Results and Discussions 

The extended size and duration of these experiments make it possible to 
follow the two-stream instability through several stages of development; first 
exponential growth, then non-linear saturation, followed by pattern growth 
and downstream advection. We identify the mechanisms responsible for these 
stages below. 

3.3.1 Magnetic Field Generation, Pattern Growth 
and Field Transport 

Encountering the shock front the incoming electrons are rapidly (being lighter 
than the ions) defie cted bv field fiuctuat i ons gr owing as a result of the two- 



grow non-linear as the defiected electrons collect into first caustic surfaces 
and then current channels fFig. l3.l() . Both streaming and rest frame electrons 
are defiected, by arguments of symmetry. 

In accordance with Ampere's law the current channels are surrounded 
by approximately cylindrical magnetic fields (illustrated by arrows in Fig. 
13.11) . causing mutual attraction between the current channels. The current 
channels thus merge in a race where larger electron channels consume smaller, 
neighbouring channels. In this manner, the transverse magnetic field grows 
in strength and scale downstream. This continues until the fields grow strong 
enough to defiect the much heavier ions into the magnetic voids between the 
electron channels. The ion channels are then subjected to the same growth 
mechanism as the electrons. When ion channels grow sufficiently powerful, 
they begin to experience Debye shielding by the electrons, which by then have 
been significantly heated by scattering on the growing electromagnetic field 
structures. The two electron populations, initially separated in 7v-space, 
merge to a single population in approximately 205e {z = 80-200) as seen in 
Fig. 13.61 The same trend is seen for the ions - albeit at a rate slower in 



stream 




The initial perturbations 
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Figure 3.2: Electron (top) and ion (bottom) currents, averaged over the y- 
direction, at time t = 1200 (240a;~/). However, due to the finite propagation time 
of the particles, these two figures can be interpreted as earlier and later times. 



proportion to mj/mg. 

The Debye shielding quenches the electron channels, while at the same 
time supporting the ion-channels; the large random velocities of the electron 
population allow the concentrated ion channels to keep sustaining strong 
magnetic fields. Figure 13.11 shows the highly concentrated ion currents, the 
more diffuse - and shielding - electron currents, and the resulting magnetic 
field. The electron and ion channels are further illustrated in Fig. 13.21 Note 
the limited z-extent of the electron current channels, while the ion current 
channels extend throughout the length of the box, merging to form larger 
scales downstream. Because of the longitudinal current channels the mag- 
netic field is predominantly transversal; we find |i?z|/|i?tot| ~ 10~^ — 10~^. 

Figure 13.31 shows the temporal development of the transverse magnetic 
field scales around z = 250. The power spectra follow power-laws, with the 
largest scales growing with time. The dominant scales at these z are of the 
order 6i at early times. Later they become comparable to L^ y. Figure EiU 
captures this scaling behavior as a function of depth for t = 2400 (480ci;~g^). 

The time evolutions of the electric and magnetic field energies are shown 
in Fig. 13.51 Seeded by fiuctuations in the fields, mass and charge den- 
sity, the two-stream instability initially grows super-linearly (t = 80 — 100 
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Figure 3.3: Power spectrum of B_|_ for z = 250 at different times. 




z 



Figure 3.4: Relative electromagnetic energy density e^. The contour color plot 
shows the power in the transverse magnetic field through the box distributed on 
spatial Fourier modes at t = 2400 (480u;~g^), with the dotted line marking the 
wavenumber with maximum power. Superposed is the spatial distribution of e^, 
averaged across the beam, at t = 2320 (464a;-e^) (dashed-dotted) and t = 2400 
(240u;~g^)(full drawn), highlighting how EM-fields are advected down through the 
box. 
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Figure 3.5: Total magnetic (full drawn) and electric (dashed) energy in the box 
as a function of time. The inset shows a log-log plot of the same data. 



= 16 — 20lu~^), reflecting approximate exponential growth in a small sub- 
volume. Subsequently the total magnetic energy grows more linearly, re- 
flecting essentially the increasing volume filling factor as the non-linearly 
saturated magnetic field structures are advected downstream. 

At t ~ 1100 the slope drops off, because of advection of the generated 
fields out of the box. The continued slow growth, for t > 1100 (220u;~g^), 
reflects the increase of the pattern size with time (cf. Fig. 13. 3|) . A larger 
pattern size corresponds to, on the average, a larger mean magnetic energy, 
since the total electric current is split up into fewer but stronger ion current 
channels. The magnetic energy scales with the square of the electric current, 
which in turn grows in inverse proportion to the number of current channels. 
The net effect is that the mean magnetic energy increases accordingly. 

The magnetic energy density keeps growing throughout our experiment, 
even though the duration of the experiment (480 u}pe) significantly exceeds 
the particle crossing time, and also exceeds the advection time of the mag- 
netic field s tructures through the box. This is in contrast to the results 



reported by ISilva et al 



( 2003[ ). where the magnetic energy density decays 
It is indeed obvious from the preceding discussion 
that the ion-electron asymmetry is essential for the survival of the current 



after about 10-30 uj„J^. 
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channels. 

From the requirement that the total plasma momentum should be con- 
served, the (electro)magnetic field produced by the two-stream instability 
acquires part of the z-momentum lost by the two-stream population in the 
shock; this opens the possibility that magnetic field structures created in the 
shock migrate downstream of the shock and thus carry away some of the 
momentum impinging on the shock. 

Our experiments show that this does indeed happen; the continuous in- 
jection of momentum transports the generated field structures downstream 
at an accelerated advection speed. The dragging of field structures through 
the dense plasma acts as to transfer momentum between the in-streaming 
and the shocked plasmas. 



3.3.2 Thermalization and Plasma Heating 

At late times the entering electrons are effectively scattered and thermalized: 
The magnetic field isotropizes the velocity distribution whereas the electric 
field generated by the e~-p charge separation acts to thermalize the popula- 
tions. Figure Eini shows that this happens over the ~ 20 electron skin depths 
from around 2; = 80 - 200. The ions are expected to also thermalize, given 
sufficient space and time. This fact leaves the massive ion bulk momentum 
constituting a vast energy reservoir for further electron heating and acceler- 
ation. Also seen in Fig. 13.61 the ions beams stay clearly separated in phase 
space, and are only slowly broadened (and heated). 

We do not see indications of a super-thermal tail in the heated electron 
distributions, and there is thus no sign of sec ond order Fermi- a cceler ation 
in the exp e rimeii t presented in this chapter. Nishikawa et al. ( 2003[ ) and 



Silva et al.l (l2003t ) reported acceleration of particles in experiments similar 



to the current experiment, except for more li mited sizes and durations, and 
the use of an e~e~^ plasma ( Silva et al. . 2003| ). On closer examination of the 



published results it appears that there is no a ctual disagreeme n t rega rding 
the absence of accelerated particles. Whence, iNishikawa et al. ( 20031 ) refer 



to transversal velocities of the order of 0.2c (their Fig. 3b), at a time where 
our experiment shows similar transv ersal velocit i es (cf. Fig. 13. 6|) that later 
develop a purely thermal spectrum. ISilva et aD (j2003[ ) refer to transversal 
velocity amplitudes up to about 0.8c (their Fig. 4), or vy ~ 2, with a shape 
of the distribution function that appears to be compatible with thermal. 
In comparison, the electron distribution illustrated by the scatter plot in 
Fig. 13.61 covers a similar interval of vj, with distribution functions that are 
close to (Lorentz-boosted) relativistic Maxwellians. Thus, there is so far 
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Figure 3.6: Thermalization and longitudinal acceleration, illustrated by scatter 
plots of the electron (orange) and ion (blue) populations. Note the back-scattered 
electron population (vz^iv) < 0). 



no compelling evidence for non-thermal particle acceleration in experiments 
with no imposed external magnetic field. Thermalization is a more likely 
ca use of the increases in tra nsversal velocities. 

Frederiksen et al. ( 2nn2t ) reported evidence for particle acceleration, with 



electron gamma up to ~ 100, in experiments with an external magnetic 
field present in the up-stream plasma. This is indeed a more promising 
sc enario for particl e acce leration experiments (although in the experiments 
bv lNishikawa et alLliooi results with an external magnetic field were similar 
to those without). Figure shows the presence of a population of back- 
scattered electrons {vz^ < 0). In the presence of an external magnetic field 
in the in-streaming plasma, this possibly facilitates Fermi acceleration in the 
shock. 



3.4 Conclusions 



The experiment reported upon here illustrates a number of fundamental 
properties of relativistic, coUisionless shocks: 
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1. Even in the absence of a magnetic field in the up-stream plasma, 
a small scale, fluctuating, and predominantly transversal magnetic field is 
unavoidably generated by a two-stream instability reminiscent of the Weibel- 
inst ability. In the current experiment the magnetic energy density reaches a 
few percent of the energy density of the in-coming beam. 

2. In the case of an e~p plasma the electrons are rapidly thermalized, 
while the ions form current channels that are the sources of deeply pene- 
trating magnetic field structures. The channels merge in the downstream 
direction, with a corresponding increase of the average magnetic energy with 
shock depth. This is expected to continue as long as a surplus of bulk relative 
momentum remains in the counter-streaming plasmas. 

3. The generated magnetic field patterns are advected downstream at 
speeds intermediate between the in-coming plasma and the rest frame plasma. 
The electromagnetic field structures thus provide scattering centers that in- 
teract with both the fast, in-coming plasma, and with the plasma that is 
initially at rest. As a result the electron populations of both components 
quickly thermalize and form a single, Lorentz-boosted thermal electron popu- 
lation. The two ion populations merge much more slowly, with only gradually 
increasing ion temperatures. 

4. The observed strong turbulence in the field structures at the shocked 
streaming interface provides a promising environment for particle accelera- 
tion. 

We emphasize that quantification of the interdependence and develop- 
ment of eu and is accessible by means of such experiments as reported 
upon here. 

Rather than devising abstract scalar parameters and eu-, that may be 
expected to depend on shock depth, media densities etc., a better approach is 
to compute synthetic radiation spectra directly from the models (see Chapter 
[7j), and then apply scaling laws to predict what would be observed from 
corresponding, real supernova remnants and Gamma-Ray Burst afterglow 
shocks. 



Chapter 4 



Large-scale Two-Dimensional 
Simulations 

4.1 Introduction 

Large-scale simulations of coUisionless electron-proton plasma shocks that are 
covering the full shock ramp are crucial for our interpretation of the radiation 
that we receive from relativistic jets (e.g. gamma-ray bursts). It is important 
to understand the nature and interdependencies of the shock "parameters" 
ee and eb- The complexity and non-linearity of coUisionless shocks give PIC 
code simulations a central role in the process of seeking further progress in 
our understanding of these shocks. 

Even though the continuously increasing computational power has reached 
a level where we may now begin to explore the full, three-dimensional col- 
lisionless shock problem, it is still only barely possible. The 3D simulations 
that were presented in Chapter |21 take from weeks to months to run on a 
modern supercomputer (parallized using 8-32 processors and 100 Gb of mem- 
ory). Still it is not possible to resolve the full shock transitions region for 
an electron-proton coUisionless plasma shock. Some of the great unanswered 
questions that cannot be targeted yet, because of this limitation in computer 
power, are for example: How large do the structures that are created by 
merging of current filaments grow in the non-linear phase of the Weibel two- 
stream instability? How many ion skin depths is the shock transition region? 
What is the fate of the generated magnetic structures in the downstream re- 
gion where the ions are thermalized and the instability has saturated? What 
is the ion and electron distribution function behind the shock ramp? Does 
the Fermi acceleration mechanism work in coUisionless shocks? Etc. 



43 



44 CHAPTER 4. LARGE-SCALE TWO-DIMENSIONAL SIMULATIONS 



4.2 Simulations and results 

To probe for answers to some of these questions, it is tempting to return to 
only two spatial dimensions. To test whether 2D simulations can account for 
the Weibel two-stream instability, we have made a 2D3V (two spatial and 
three velocity dimensions - sometimes referred to as 2|D) repetition of the 
simulation described in Chapter El We use the exact same plasma densities 
and velocities as in the 3D experiment: In the shock rest frame, ISM particles 
are injected with F = 3 into a plasma, initially at rest, with electron plasma 
frequency Upe = 0.2 c/A (simulation units). The major difference is that 
the simulation box is two-dimensional with 800 x 4000 grid zones instead of 
the 3D 200 x 200 x 800. The 2D box size corresponds to 40 x 200 ion skin 
depths. Compared to the 3D box, the extra grid zones in the ^-direction is 
used to cover 320 grid zones upstream of the 3D box and and 2880 grid zones 
downstream of the 3D box. 
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Figure 4.1: A comparison of density {top) and magnetic field generation (bottom) 
between 2D (dashed blue) and 3D (solid red). The three columns corresponds to 
snapshots taken at time t = 40, 80 and 120 <j-'~/. 

Figure \AA\ shows the evolution of the plasma density (top panel) and mag- 
netic energy (bottom panel) in the 2D and 3D simulation. A good agreement 
between the two simulations is found, at least with respect to the evolu- 
tion of the density. The minor differences are mainly caused by the artificial 
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reflection of thermal particles at the leftmost boundary, and the drain of par- 
ticles at the rightmost boundary in the 3D simulations. For the generation of 
magnetic field, the situation is somewhat different. In the very early stage, 
the field generated in the 3D simulation is stronger than in 2D, because of 
the extra transverse dimension that the 3D instability can collect particles 
from. However, at later times the growth in magnetic field in the 2D run 
exceeds the 3D case. This is primarily caused by two effects: 1) in the 2D 
box, upscattered particles (defined as shocked particles that are traveling 
upstream into the ISM) can generate a seed field further upstream. Such 
a seed field enhances the growth of the instability. This mechanism cannot 
be followed (yet) in the smaller 3D box. 2) ion current channels can merge 
to larger transverse structures. These also cannot be followed in the smaller 
3D box. To a first approximation, the two-dimensional simulation shows a 
promising resemblance with the three-dimensional simulations, and poten- 
tially can surpass them in some respects. Encouraged by these results, we 
have continued the large-scale 2D simulations in order to follow the long time 
evolution (tmax = 400u;~/). This duration corresponds to two light crossing 
times of this larger box. 




z (ion skin depth) 

Figure 4.2: The density profile at different time, sampled in the time range 
t = — 4:00uj~-^ with an interval of 20uj~-^ . The different time snapshots are color 
coded with red for early times through blue to black at late times. 



The evolution of the plasma density as function of shock depth is shown 
in Fig. 14.21 for different times. The time-range is t = - 400u;"/ with an 
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interval of 200;^^ • is immediately clear that the initial setup is not con- 
sistent with a steady state solution. This is not surprising, since initially 
there are no electromagnetic fields to facilitate the momentum transfer be- 
tween the injected and quiescent material. As the simulation evolves, the 
shock converges toward the real physical solution. In three-dimensional sim- 
ulations performed by us and other groups, it has been neither temporally 
nor spatially possible to follow the evolution to this stage. However, the 2D 
simulations that we have performed cover a much larger range in space and 
time. 

At late times, the shock reaches a steady state where the value of the main 
density peak is constant in time. The whole shock structure is propagating 
with 0.35c in the downstream direction. Such a steady state is expected 
when the computational domain covers the full shock transition ramp. In 
this case, the injected ISM plasma should have thermalized and merged with 
the shock plasma. Thus, only a single plasma population should be present 
in the downstream region. This is indeed the case, as seen in Fig. 14.31 The 
figure shows the ion [left) and electron (right) distribution functions of the 
momentum along the shock propagation direction (7f^), at time t = 400ci;~g^. 
The four different curves correspond to different regions of the computational 
domain ( ^ = [0, 50] , ^ = [50,100], z = [100,150], and z = [150,200] ion 
skin depths). In the shock ramp [z = [0, 150] ion skin depths), where ther- 
malization of the two ion populations (ISM and shocked) has not occurred 
yet, two distinct peaks are seen. In the last segment {z = [150, 200] ion skin 
depths), there exist only one single peak. This means that the ions are fully 
thermalized. It is also noteworthy that the merged ion population translates 
downstream with '-fVz — 0.4, consistent with the estimate found from figure 
14.21 Even though the peak moves relatively slowly, the particles that define 
the peak are moving relativistically and are continuously getting replaced 
with "fresh" particles. 

Two shock parameters of great interest are the equipartition parameters 
ee and e^- They describe the amount of energy that is deposited in the 
electrons and in the magnetic field, relative to the kinetic energy of the 
ions. These parameters are normally treated as free parameters, available to 
cover our general ignorance of the details of the microphysics in coUisionless 
shocks. Fig. 14.41 shows eg (left) and es (right) as functions of the shock depth. 
It is found that the electrons reach close to equipartition with the ions. 
This, combined with the arguments above, allows us to conclude that the 
simulations cover the full shock. In the downstream region, equipartition has 
not yet been reached. This is because since everything moves with velocities 
close to the speed of light, we see earlier stages in the evolution the further 
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Figure 4.3: The ion (left) and electron {right) distribution functions of the 
momentum {■jVz) along the shock propagation-direction. The particles are sampled 
at t = 4000;^^"'^ in four different segments of the shock ramp: z = [0,50] {blue, 
dotted), z = [50,100] {black, dashed), z = [100,150] {green, dot-dashed) and z = 
[150,200] {red, solid). The (thin black dotted) vertical lines indicate = 
and the injection momentum (corresponding to 7 = 3). Note that in the last 
segment {z = [150, 200] ion skin depths) the ions are very close to being one single 
population. This indicates that the full shock ramp is included in the box. 



downstream we look. The 63 plot shows us that most of the field generation 
takes place in the foreshock. This is not surprising, since the streaming 
instability in its nature requires counter-streaming populations (anisotropy). 
Downstream, the plasma populations are fully thermalized into on single 
population. 

Finally, to give an impression of the extremely complex and non-linear na- 
ture of the coUisionless shock, we present contour plots of the large-scale evo- 
lution of the shock in Figure 1^31 The figure shows the ion density structures 
at four different time-steps {t = 100,150,200, and 300u;~/). At early times 
we see that the Weibel two-stream instability is occurring semi-symmetrically 
at the contact discontinuities both upstream and downstream. 



4.3 Summary 

Even with the computational power available on supercomputers today, it 
is not possible to make three-dimensional PIC simulations that cover a full 
coUisionless shock, while also resolving the microphysics. This has lead us 
to investigate whether two-dimensional simulations can include the physical 
processes that dominate these shocks (mainly the non-linear Weibel instabil- 
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Figure 4.4: Electron (left) and magnetic {right) equipartition values (eg and es) 
at time t = SOOo;"-^. The dotted line in both plots shows the shape of the total 
plasma density for comparison. The electrons reach up to 80% of equipartition 
and the magnetic field up to 10%. 



ity). We find that this is the case as far as 3D results exist to compare with. 
The density profile and generation of magnetic fields are quite similar in the 
two cases. Encouraged by this, we have made large-scale 2D PIC simulations 
of the formation of a collisionless shock, for an ejecta propagating with a bulk 
lorentz factor P = 3. The simulation box cove rs 40 x 200 ion skin de pths 
(compared to the 10 x 10 x 4 in the 3D case bvlFrederiksen et al.ll2004L and 
4 X 4 X 15 in simulations by Nishikawa et al. 2005[ l. 

In the results we find that the 2D simulations cover a spatial region just 
large enough for the ISM plasma (injected into a plasma in the shock rest- 
frame) to thermalize and merge with the shocked plasma. Even though the 
initial conditions are not setup correctly, the system converges to a quasi- 
stationary state where at least three different regions are identified: 1) An 
upstream foreshock with great anisotropy, consisting of in-streaming ISM 
plasma and upscattered shock particles. In this region, strong magnetic field 
generation is taking place with = 5 — 10%. 2) A dense, hot thermalization 
region where the electrons and ions are close to equipartition (eg = 50 — 80%). 
In this region, the magnetic field is relatively weak (e^ ^ 1%) but still strong 
enough to account for most GRB afterglow estimates. 3) A hot downstream 
region where the magnetic field is of the order percents of equipartition. This 
region is clumpy, with structure sizes of the order 1 — 20 ion skin depths (see 
Fig. 1131). 

In these 2D simulations it is found that with an initial condition where 
a dense plasma is at rest in the simulation frame, the resulting shock rest 
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Figure 4.5: Plasma density contour plots at four different time levels A, B, C 
and D corresponding to t = 100, 150, 200 and 300 uj~l 



frame was propagating with 0.35c in the downstream direction. Based on 
these results, we plan to rerun the simulations in a frame where the dense 
population initially travels with 0.35c in the upstream direction. This would 
potentially allow us to follow the shock for longer periods of time. 

Finally we note that the upscattering of particles might be important for 
Fermi acceleration. It is likely that these particles can be deflected back into 
the shock region. To test whether Fermi acceleration is taking place requires 
even larger simulations and the question can only be answered by future 
work. 



Chapter 5 

Shocks in Ambient Magnetic 
Fields 

We now continue with an investigation of the Weibel i nstabihty in coUision- 



less s hocks in the presence of an ambient magnetic field (jHededal and Nishikawa 



2nn5[ l. This scenario is important especially in the internal shock phase and 



in the very early afterglow, where it is possible that a magnetic field, carried 
from the GRB progenitor, still exist in the ejecta. 

Plasma outflows from gamma-ray bursts, supernovae, and relativistic jets, 
in general interact with the surrounding medium through coUisionless shocks. 
The microphysical details of such shocks are still poorly understood, which, 
potentially, can introduce uncertainties in the interpretation of observations. 
It is now well established that the Weibel two-stream instability is capable 
of generating strong electromagnetic fields in the transition region between 
the jet and the ambient plasma. However, the parameter space of coUi- 
sionless shocks is vast and still remains unexplored. In this chapter, we 
focus on how an ambient magnetic field affects the evolution of the elec- 
tron Weibel instability and the associated shock. Using a particle-in-cell 
code, we have performed three-dimensional numerical experiments on such 
shocks. We compare simulations in which a jet is injected into an unmag- 
netized plasma with simulations in which the jet is injected into a plasma 
with an ambient magnetic field both parallel and perpendicular to the jet 
flow. We find that there exists a threshold of the magnetic field strengths, 
below which the Weibel two-stream instability dominates, and we note that 
the interstellar medium magnetic field strength lies well below this value. In 
the case of a strong magnetic field parallel to the jet, the Weibel instability 
is quenched. In the strong perpendicular case, ambient and jet electrons are 
strongly accelerated because of the charge separation between deflected jet 
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electrons and less deflected jet ions. Also, the electromagnetic topologies 
become highly non-linear and complex with the appearance of anti-parallel 
fleld conflgurations. 



5.1 Introduction 



The coUisionless plasma condition applies to many astrophysical scenarios, 
including, the outflow from gamma-ray bursts (GRBs), active galactic nuclei, 
and relativistic jets in general. The complexity of kinetic effects and instabil- 
ities makes it difficult to understand the nature of coUisionless shocks. Only 
recently, the increase in available computational power has made it possible 
to investigate the full three-dimensional dynamics of coUisionless shocks. 

In the context of GRB afterglows, observations indicate that shock-compressed 
magnetic field from the interstellar medium (ISM) is several orders of magni- 
tude too weak to match observations. Particle-in-cell (PIC) simulations have 
revealed that the Weibel two-stream instability is capable of generating the 
required el ectromagnetic field strength of the order of percents of equiparti- 
tion value dKazimura et al" ^ 19981: Medvedev and Loeb ^ 19991: Nishikawa et al. , 



20031 : iNishikawa et all 1200,4 ISilva et alLbnoilFrederiksen et al.L 120041 ). Fu7 



thermore, PIC simulations have shown that in situ n on thermal particle ac- 
celeration takes place in the shock transition regio n (jHoshino and Shimada . 
2OO2I: iHededal et al.1 . |2004 [ ISaito and Sakail . l2003 [) . Three-dimension a l sim- 
ulations using ~ 10^ electron-positron pairs by ISakai and Matsuol (l2004l ) 



showed how complex magnetic topologies are formed when injecting a mildly 
relativistic jet into a force-free magnetic field with both p arallel and per- 



pendic ular components. With a two-dimensional analysis, ISaito and Sakai 
( 20031 ) lound that an ambient parallel magnetic field can quench the two- 
stream instability in the weakly relativistic case. In this chapter, we use 
three-dimensional PIC experiments to investigate how the two-stream insta- 
bility is affected by the presence of an ambient magnetic field. Using up 
to ~ 10^ particles and 125 x 125 x 1200 grid zones, we investigate the de- 
velopment of complex magnetic topologies when injecting a fully relativistic 
jet (bulk Lorentz factor P = 5) into an ambient magnetized plasma. Using 
varying field strengths, we focus on the case of a transverse magnetic field 
and compare it with the case of a parallel magnetic field. The simulations 
are mainly concerned with the electron dynamics since processes involving 
the heavier ions evolve on much longer timescales. 



5.2. THE NUMERICAL EXPERIMENTS 



53 



5.2 The numerical experiments 



We use the PIC code described bv iFrederiksen et al.l ()2004 ). The code works 



from first principles and evolves the equation of motion for the particles 
together with Maxwell's equations. In the simulation experiments, we inject 
an electron-proton plasma (a jet) into an ambient plasma (the ISM) initially 
at rest (Fig. 15. The jet is moving with a relativistic velocity of 0.98c 
along the z-direction, corresponding to Lorentz factor '-jjet = 5. The ion-to- 
electron mass ratio is set to mi/mc = 20. The jet plasma and the ambient 
plasma have the same density, n, and the corresponding electron plasma rest- 
frame frequency ojpe = [ne^/{'meeo)] = 0.035A^^ (e is the unit charge, eo the 
vacuum permittivity, and At the simulation time unit). We choose this low 
value in order to properly resolve the microphysics. Initially, the interface 
between the ambient and the injected plasma is located at 2; = 3Ae, where Ae 
is the electron skin depth defined as Ae = c/upe = 28. 6 (c is the speed of 
light, and A^; the grid size). The time step and grid size obey the Courant- 
Friedrichs-Levy condition A^ = 0.5Ax/c. Both plasma populations are, in 
their respective rest frames, Maxwellian distributed with a thermal electron 
velocity Vth — 0.03c. This temperature allows us to numerically resolve the 
plasma Debye length with at least one grid length. 




Figure 5.1: Schematic example of the simulation setup. Jet plasma is homoge- 
neously and continuously injected in the z-direction throughout the x — y plane 
at z = 0. Inside, the box is populated by a plasma population, initially at rest. 
In this specific example, an ambient magnetic field is set up in the x-direction 
(perpendicular case). 



We consider three different ambient magnetic configurations: no magnetic 
field, a magnetic field parallel to the flow, and a magnetic field perpendicular 
to the flow. The magnetic fleld is initially setup to be homogeneous and at 
rest in the ambient plasma. The experiments are carried out with 1 billion 
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particles inside 125 x 125 x 1200 grid zones. In terms of electron skin depths, 
this corresponds to 4.4 x 4.4 x 42Ae. The boundary conditions are periodic 
in the direction transverse to the jet flow {x, y). In the parallel direction, jet 
particles are continuously injected at the leftmost boundary (2 = 0). At the 
leftmost and rightmost z boundary, electromagnetic waves are absorbed, and 
we allow particles to escape in order to avoid unphysical feedback. The total 
energy throughout the simulations is conserved with an error less than 1%. 



5.3 Results 



Initially, we ran simulations with no ambient magnetic field and observed 
the growth of the Weib e l two-stream instabil it y also found in previou s work 



dKazimura et al" . 19981: Medvedev and Loeb . 19991: Nishi cawa et all . I2OO3I : 



Nishikawa et al. ri200,4ISilva et all . 1200^1: IPrederiksen et al.L 



20041 ). The Weibel 



two-stream instability works when magnetic perturbations transverse to the 
flow collect streaming particles into current bundles that in turn amplify the 
magnetic perturbations. In the non-linear stage, we observe how current fil- 
aments merge into increasingly larger patterns. The electromagnetic energy 
grows to €3 — 1%, where describes the amount of total injected kinetic 
energy that is converted to magnetic energy. 



5.3.1 Parallel Magnetic Field 

In the presence of a strong magnetic field component parallel to the flow, 
particles are not able to collect into bundles, since transverse velocity com- 
ponents are deflected. We have performed five runs, with parallel magnetic 
fields corresponding, respectively, to ujpf,/ujce =40, 20, 10, 5, and 1, while 
keeping tUpe constant; ujce = / {'jjetme) is the jet electron gyrofrequency. 
The resulting field generation efficiency can be seen in Fig. 15.21 at t = 21u!~^ 
where the jet front has reached z = 23Ae. In the case of ujpe/oJce = 40, the 
Weibel instability overcomes the parallel field, and although initially slightly 
suppressed, it eventually evolves as in the case of no ambient magnetic field. 
Increasing the magnetic field to Upe/uce = 1 effectively suppresses the insta- 
bihty. Thus, for an ISM strength magnetic field {ujpe/uce — 1500) parallel to 
the plasma flow, the Weibel two-stream instability evolves unhindered, and 
the generated field will exceed the ISM field. We find from the simulations 
that it would take a milligauss strength parallel magnetic field to effectively 
quench the instability for a 7 = 5 jet expanding in an environment with 
density similar to the ISM. 
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z (electron skindepth) 

Figure 5.2: Growth of the Weibel two-stream instabihty for different strengths 
of the parallel ambient magnetic field at time t = 21lo~^. Here we show the effec- 
tiveness of the field generation, as measured by the average transverse magnetic 
field amplitude as a function of z. The solid line corresponds to ujpe/uJce = 40, 
the dotted line to 20, the dashed line to 10, the dot-dashed line to 5, and the 
triple-dot-dashed line to 1. The case of no ambient magnetic field is very similar 
to that of LOpe/uJce = 40. The magnetic field amplitude is in arbitrary units 



The left panel of Fig. 15.31 shows the resulting electron momentum dis- 
tribution function for different values of Upe/uJce at t = 30u!~^. Since the 
presence of a strong parallel magnetic field suppresses the generation of a 
transverse magnetic field, there exists no mechanism that can heat the elec- 
trons and transfer momentum between the two electron populations. Thus 
the jet plasma propagates unperturbed. Where there is no parallel magnetic 
field or only a weak magnetic field {upe/uce = 1500), we observe how the jet 
and ambient plasma is heated and how momentum is transferred between 
the two populations. 

5.3.2 Perpendicular Magnetic Field 

We have performed experiments with an ambient magnetic field perpendicu- 
lar to the jet flow (Fig.EH]), with field strengths corresponding to ujpe/ujce=^500, 
40, 20 and 5. By including the displacement current, one can derive the rela- 
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Figure 5.3: Normalized electron momentum distribution functions at time 30uj~^. 
The left panel is for runs with an ambient magnetic field parallel to the injected 
plasma with Wpe/wce = 1 {solid line), 20 {dotted line), 40 {dashed line), and 1500 
{dot-dashed line). The right panel is for runs in which the initial magnetic field is 
perpendicular to the inflow and ujpe/uJce = 5 {solid line), 20 {dotted line), 40 {dashed 
line) and 1500 {dot-dashed line). The vertical line shows the injected momentum 
7 = 5. The distribution functions are for electrons with z > 15Ae. 



tivistic Alfven speed t;;^^ _ c-2+(t;7"-'^'='-)-2 ^ c/[l+(^peMe)^K/me)77J]^/2^ 
where y'^^-'^'^'-- = B/[fj.Qn{mi + m^^l'^ is the non-relativistic counterpart. 
From this we calculate the corresponding relativistic Alfven Mach numbers 
IjetV.etlvA = 6572, 175, 88, and 22. 

Again, the Upe/uce = 1500 run has been chosen because it resembles 
the typical density and microgauss magnetic field strength of the ISM. We 
find that the magnetic field generated by the two-stream instability domi- 
nates the ambient magnetic field, and the result resembles the case with no 
ambient magnetic field. Furthermore, as expected in both the parallel and 
perpendicular cases, the electron momentum distributions (Fig. 15. 3|) are very 
similar, except for a weak merging between the ambient and jet electrons in 
the perpendicular case. 

In the run with uOpe/oJce = 20, the result differs substantially from the 
previous cases. With reference to Fig. 15.41 we describe the different stages 
of the evolution: Initially, the injected particles are deflected by the ambient 
magnetic field. The magnetic field is piled up behind the jet front, and the 
enhanced magnetic fields bend jet electron trajectories sharply. This has two 
implications. First, the ions, being more massive, penetrate deeper than the 
deflected electrons. This creates a charge separation near the jet head that 
effectively accelerates both ambient and injected electrons behind the ion jet 
front as shown in Fig. 15.51 Second, the deflected electrons eventually become 
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Figure 5.4: Snapshot at t = 16ti;~g of the highly comphcated magnetic field 
topology resulting when a jet is injected into a plasma with an ambient magnetic 
field transverse to the jet flow. The bottom panel shows magnetic field lines in a 
subsection of the computational box (from z = 9Ae to z = ISAg. The top panel 
refers to a schematic explanation in the x — z plane: Jet electrons are bent by the 
ambient magnetic fleld (region A). As a result of the Weibel instability, the elec- 
trons bundle into current beams (region C) that in turn reverse the field topology 
(region B). This will eventually bend the jet beam in the opposite direction. 



subject to the Weibel two-stream instability. This forms electron current 
channels at some angle to the initial direction of injection, as shown in the 
upper panel of Fig. 15.41 Around these current channels, magnetic loops are 
induced (Fig. 15.41 region C). Magnetic islands are formed and the ambient 
magnetic field is reversed behind the loops (region B). In this region, we 
find acceleration of electrons in the x-direction. The activity in this region 
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has similarities to reconnection, but it is beyond the scope of this thesis to 
investigate this topic. 

Electrons Ions 
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Figure 5.5: Phase-space plot of electrons {left column) and ions {right column) at 
time t = 30uj~} . The data is from a simulation of a shock with transverse ambient 
magnetic field corresponding to ujpe/uJce = 20. The three rows correspond to 'jVx, 
jVy and ^Vz- The jet plasma (black dots) are injected at z = with 7 = 5. The 
ambient electrons {red dots) are initially at rest but are strongly accelerated by 
the jet. 



In other regions, the ambient magnetic field is strongly compressed, and 
this amplifies the field strength up to 5 times the initial value. As a result. 
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parts of the jet electrons are reversed in their direction. This may be seen 
from Fig. 15.51 which shows a phase-space plot of both ambient and jet elec- 
trons at t = 30a;~g^. We see several interesting features here. In the region 
2; = (15 — 21) Ae, we observe how ambient electrons are swept up by the jet. 
Behind the jet front, both ambient and jet electrons are strongly accelerated 
since the jet ions, being heavier, take a straighter path than the jet electrons, 
and this creates a strong charge separation. The excess of positive charge at 
the very front of the jet head is very persistent and hard to shield since the 
jet ions are moving close to the speed of light. Thus, there is a continuous 
transfer of z-momentum from the jet ions to the electrons. In the case of the 
perpendicular ambient field, more violent processes take place than in the 
case of the parallel field, which may be seen in Fig. 15. HI [right panel). Here 
we see that mixing of the two plasma populations is much more effective 
for the perpendicular case. However, the spectrum of the electron's momen- 
tum is highly nonthermal, with strong acceleration of both jet and ambient 
electrons. The cutoff in electron acceleration depends on the magnetic field 
strength. The maximum {'^v\\ ~ 10) aX z = 20Ae in Fig 15.51 corresponds to 
the cutoff shown by the dotted line in the right panel in Fig. 15.31 It should 
be noted that the current channels that are caused by the bent jet electron 
trajectories at the early time, as shown in Fig. 15.41 are also seen in Fig. 15.51 
The first current channels have moved to around z = 20Ae. At z = 15Ae, a 
second current channel is created by the deflected jet electrons. This periodic 
phenomenon involves the ions (Fig. 15. 6p in a highly nonlinear process but is 
beyond the scope of this thesis and will be investigated in subsequent work. 



5.4 Conclusions 



Using a three-dimensional relativistic particle-in-cell code, we have investi- 
gated how an ambient magnetic field affects the dynamics of a relativistic jet 
in the coUisionless shock region. We have examined how the different ambient 
magnetic topology and strength affect the growth of the electron Weibel two- 
stream instability and the associated electron acceleration. This instability 
is an important mechanism in coUisionless shocks. It facilitates rnomen- 

tum transfer between coll i ding plasma populations dKazimura et al. j_ 19981: 

Medvedev and LoebLll999l:lNishikawa et ahlEinilNishikawa et al.l . 1200,4 ' Silva et al 
2OO3I : Frederiksen et al.L 12004 ) and can accelerate electrons to nonthermal dis- 
tribu tions ( Hoshino and Shimada , 20021 : Hededal et al. , 2004 ; Saito and Sakai , 
2003[ ). CoUisionless shocks are found in the interface between relativistic 
outflows (e.g., from gamma-ray bursts and active galactic nuclei) and the 
surrounding medium (e.g., the ISM). 
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Figure 5.6: The ion density averaged over the box width as a function of shock 
depth. The different colors correspond to different times. The simulation is for 
uipe/oJce = 20. The distance between each peak is comparable to the ion skin 
depth. 



We find substantial differences between the cases of ambient magnetic 
fields transverse and parallel to the jet flow. However, common to both cases 
is that it takes an ambient magnetic field strength much stronger than the 
strength of the magnetic field typically found in the ISM to effectively sup- 
press the Weibel two-stream instability. In the case of a parallel magnetic 
field, ujpe/oJce must be smaller than 5 to effectively suppress the instabil- 
ity. This result is in g ood agreement with two-dimensional simulations by 



Saito and Sakail (j2003[ ). and thus this limit seems independent of 7jet- For 
a typical ISM density of 10^m~^, this corresponds to a milligauss magnetic 
field. We emphasize the role of lOy^jiOce as an important parameter for c oUi- 
sionless shocks, as was also pointed out by Shimada and Hoshinol (2004). 

In the case of perpendicular injection, the dynamics are different from 
the parallel injection. Here, the electrons are deflected by the magnetic 
field, and this creates a charge separation from the less deflected ions. The 
charge separation drags the ambient and jet electrons, and consequently they 
are strongly accelerated along the ^-direction. Furthermore, as a result of 
the Weibel instability, current channels are generated around the ambient 
magnetic field, which is curled and locally reversed. 

These simulations provide insights into the complex dynamics of rela- 
tivistic jets. Further investigations are required to understand the detailed 
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physics involved. Larger simulations (above 10^ particles) with longer boxes 
are needed to cover the instability domain of the ions, to investigate the 
full evolution of the complicated dynamics, and to resolve the whole shock 
ramp. 



Chapter 6 



Non— Fermi Acceleration in 
Plasma Shocks 



Collisionless plasma shock theory, which apphes, for example, to the after- 
glow of gamma-ray bursts, still contains key issues that are poorly under- 
stood. In this chapter, I discuss the results of charged particle dynamics 
i n a highly relativisti c collisionless shock numerically using ~ 10^ particles 



( Hededal et all 120041 1 . We find a power-law distribution of accelerated elec- 
trons, which upon detailed investigation turns out to originate from an accel- 
eration mechanism that is decidedly different from Fermi acceleration. Elec- 
trons are accelerated by strong filamentation instabilities in the shocked inter- 
penetrating plasmas and coincide spatially with the power-law-distributed 
current filamentary structures. These structures are an inevitable conse- 
quence of the now well-established Weibel-like two-stream instability that 
operates in relativistic collisionless shocks. The electrons are accelerated and 
decelerated instantaneously and locally: a scenery that differs qualitatively 
from recursive acceleration mechanisms such as Fermi acceleration. The 
slopes of the electron distribution power-laws are in concordance with the 
particle power law spectra inferred from observed afterglow synchrotron ra- 
diation in gamma-ray bursts, and the mechanism can possibly explain more 
generally the origin of nonthermal radiation from shocked interstellar and 
circumstellar regions and from relativistic jets. 



6.1 Introduction 

Given the highly relativistic conditions in the outflow from gamma-ray bursts 
(GRBs), the mean free path for particle Coulomb collisions in the afterglow 
shock is several orders of magnitude larger than the fireball itself. In explain- 
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ing the microphysical processes that work to define the shock, MHD becomes 
inadequate and coUisionless plasma shock theory stands imperative. In par- 
ticular, two key issues remain, namely, the origin and nature of the magnetic 
field in the shocked region and the mechanism by which electrons are accel- 
erated from a thermal population to a power-law distribution N{'-f)d'~f oc 7"^. 
Both ingredients are needed to explain the observed afterglow spectra (e.g. 
Kumar . 2000l : Panaitescu and Kumai . 2001 ). 



Regarding the origin of the magnetic field in the shocked region, observa- 
tions are not compatible with a compressed interstel lar magnetic field, which 
would be orders of magnitude smaller than needed (jOruzinov and Waxman , 
Il999t ). It has been suggested that a Weibel-like two-stream instabilitv can 
generate a magnetic field in the shocked region (iMedvedev and Loebl 19991 : 
Frederiksen et al.ll2002l: iNishikawa eTaDbOOi ISilva et al.ll2003t ). Computer 



experiments by Frederiksen et alT ( 2004| ) showed that the nonlinear stage of 
a two-stream instability induces a magnetic field in situ with an energy con- 
tent of a few percent of the equipartition value, consistent with that required 
by observations. 

Fermi acceleration (Fermi, 19491 ) has, so far, been widely accepted as the 
mechanism that provides the inferred electron accel eration. It has been em- 



Niemiec and Ostrowski 



ployed extensively in Monte Carlo simulations (e.g. 
(200J) and references therein), where it operates in conjunction with certain 
assumptions about the scattering of particles and the structure of the mag- 
netic field. The mechanism has, however, not been conclusi vely demonstrated 

to occ ur in ab initio particle simulations. As pointed out bv lNiemiec and Ostrowski 
( 2004^ . further significant advance in the study of relativistic shock particle 
acceleration is unlikely without u nderstanding the d e tailed microphysics of 
coUisionless shocks. Also, recentlv iBaring and Brabv (2004) found that par- 
ticle distribution functions (PDFs) inferred from GRB observations are in 
contradistinction with standard acceleration mechanisms such as diffusive 
Fermi acceleration. 

In this chapter, we study ab initio the particle dynamics in a coUisionless 
shock with bulk Lorentz factor F = 15. We find a new particle acceleration 
mechanism, which we present in Section 16.21 Detailed numerical results 
are presented and interpreted in Section 16. 3[ while Section 16.41 contains the 
conclusions. 



6.2 A new acceleration mechanism 



We have performed a series of numerical experiments where coUisionless 
shocks are created by two colliding plasma populations. These experiments 
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Figure 6.1: (A) Ray-traced electron paths (red) and current density (blue). The 
colors of the electron paths reflect their four-velocity according to the color table 
in the inset (B). The shadows are equivalent to the x and y projections of their 
paths. The ion current density is shown with blue colors according to the color 
table in the inset. The inset also shows the ion current density (blue) integrated 
along the x-axis with the spatial distribution of fast-moving electrons (red) over 
plotted. The main figure is an enlargement of the white square in the inset. 



are described in more detail below, but a common feature is that the elec- 
tron PDF has a high-energy tail that is power-law distributed. By carefully 
examining the paths of representative accelerated electrons, tracing them 
backward and forward in time, we have been able to identify the mechanism 
responsible for their acceleration. The acceleration mechanism, which as far 
as we can tell has not been discussed in the literature previously, works as 
follows: 

When two nonmagnetized collisionless plasma populations interpenetrate, 
current channels are formed through a Weibel-like two-stream instabilit'^ 



( Medvedey and Loebl 19991 : Frederiksen et al. 20021 : Nishikawa et al. 2| 



Silva et al. 20031 ). In the nonlinear stage of evolution of this instability, ion 



current channels merge into increasingly stronger patterns, while e l ectron s 



act to Debye shield these channels, as shown bv iFrederiksen et al.l (j200J). 



That work further showed that a Fourier decomposition of the transverse 
ion current fi laments exhibits power -law behavior, which has been recently 



H). 



confirmed bv lMedvedev et al. 

Figure 16.21 presents the electric and magnetic field in the vicinity of an 
ion current channel, generated by the Weibel two-stream instability. At dis- 
tances less than the Debye length, the ion current channels are surrounded 
by transverse electric fields that accelerate the electrons toward the current 
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Figure 6.2: A contour plot of the ion current strength with arrows overplotted 
to represent the magnetic field {left panel) and the electric field {right panel). 



channels. However, the magnetic fields that are induced around the current 
channels act to deflect the path of the accelerated electrons, boosting them 
instead in the direction of the ion flow. Since the forces working are caused 
by quasi-stationary fields, the acceleration is a simple consequence of poten- 
tial energy being converted into kinetic energy. Therefore, the electrons are 
decelerated again when leaving the current channel and reach their maximal 
velocities at the centers of the current channels. Hence, as illustrated by Fig. 
16.1b . the spatial distribution of the high-energy electrons is a direct match 
to the ion current channels and the properties of the accelerated electrons 
depend primarily on the local conditions in the plasma. 

One might argue that the near-potential behavior of the electrons, where 
they essentially must lose most of their energy to escape from the current 
channels, would make the mechanism uninteresting as an acceleration mech- 
anism since fast electrons cannot easily escape. However, this feature may 
instead be a major advantage, since it means that energy losses due to escape 
are small and that the electrons remain trapped long enough to have time 
to lose their energy via a combination of bremsstrahlung and synchrotron or 
jitter radiation. We observe that only a very small fraction of the electrons 
manage to escape, while still retaining most of their kinetic energy. This 
happens mainly at sudden bends or mergers of the ion channels, where the 
electron orbits cannot be described in terms of a particle moving in a static 
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electromagnetic field. 




I 



Figure 6.3: Ion current channel surrounded by an electric and a magnetic field. 
A schematic drawing of the scenario in Fig. 16.21 seen from the side. Electrons in 
the vicinity of the current channels are thus subject to a Lorentz force with both 
an electric and magnetic component, working together to accelerate the electrons 
along the ion flow. Crossing the center of the channel, the process reverses, leading 
to an oscillating movement along the channel. 



To analyze the acceleration scenario quantitatively, we use the sketch 
in Fig. 16. H[ We assume that the ion current channel has radius R, that 
the total charge inside the cylinder per unit length is pc, and that the ions 
stream with velocity u in the laboratory rest frame (see Fig. 16.31 and inset 
for definition of rest frames). Consider an electron with charge —q and mass 
m at a distance r from the center of the channel, initially having no velocity 
components perpendicular to the cylinder. If we analyze everything in the 
ion channel rest frame, the problem reduces to electrostatics and it is possible 
to analytically calculate the change in four-velocity of the electron when it 
reaches the surface of the cylinder. Since the electric force only works along 
the r-axis, the four-velocity along the 2;-axis of the electron is conserved in 
the ion channel rest frame. Hence, we can calculate both the total change 
in energy and the change in the different velocity components. Returning to 
the laboratory rest frame, we find 

^{iVz) electron = uA'Jelectron ■ (6.2) 

The change in the Lorentz boost is directly proportional to the total charge 
inside the channel and inversely proportional to the electron mass. Debye 
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shielding reduces the electric field further away from the ion channel, so the 
estimate above is only valid for distances smaller than a Debye length. 



6.3 Computer Experiments 

The experiments were performed with the three-dimensio nal relativistic ki- 



netic a nd electromagnetic particle-in-cell code described bv lFrederiksen et al 



( 20041 ) ■ The code works from first principles, by solving Maxwell's equations 
for the electromagnetic fields and solving the Lorentz force equation of mo- 
tion for the particles. 

Two colliding plasma populations are set up in the rest frame of one of the 
populations (downstream, e.g., a jet). A less dense population (upstream, 
e.g. the interstellar medium) is continuously injected at the left boundary 
with a relativistic velocity corresponding to a Lorentz factor F = 15. The 
two populations initially differ in density by a factor of 3. We use a com- 
putational box with 125 x 125 x 2000 grid points and a total of 8 x 10® 
particles. The ion rest-frame plasma frequency in the downstream medium 
is ujpi = 0.075, rendering the box 150 ion skin depths long. The electron rest- 
frame plasma frequency is ujpe = 0.3 in order to resolve also the microphysics 
of the electrons. Hence, the ion-to-electron mass ratio is rrii/me = 16. Other 
mass ratios and plasma frequencies were used in complementary experiments. 
Initially, both plasma populations are unmagnetized. 

The maximum experiment duration has t^ax = 340 cj"/, which is suf- 
ficient for the continuously injected upstream plasma (F = 15, f ~ c) 
to travel 2.3 times the length of the box. The extended size and dura- 
tion of these experiments enable observations of the streaming instabili- 
ties ar id concurrent particle a cceleration through several stages of develop- 



ment (jFrederiksen et al.Ll2004( ). Momentum losses to radiation (cooling) are 
presently not included in the model. We have, however, verified that none 
of the accelerated particles in the experiment would be subject to signifi- 
cant synchrotron cooling. The emitted radiation may thus be expected to 
accurately refiect the distribution of accelerated electrons. 

When comparing numerical data with eq. 16.11 we take r to be the ra- 
dius where Debye shielding starts to be important. Using a cross section 
approximately in the middle of Fig. 16.11 we find A^'jVz) electron = 58 In (r/i?). 
It is hard to determine exactly when Debye shielding becomes effective, but 
looking at electron paths and the profile of the electric field, we estimate that 
\n{r/R) ^ 1.3. Consequently, according to eq. 16.11 the maximally attainable 
four-velocity in this experiment is in the neighbourhood of {'yVz)max = 75. 
This is in good agreement with the results from our experiments, where the 
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Figure 6.4: A scatter plot of the local ion current density Jion versus the four 
velocity of the electrons in a region downstream of the shock. Overplotted is a line 
(thin) showing the average four velocity as a function of J/on, and a line (thick) 
showing a straight line fit. Because 'cold' trapped thermal electrons (indicated 
with the ellipse) exist inside the ion current channel they count towards lowering 
the average four velocity at high J/on- If we cleaned our scatter plot, statistically 
removing all thermal electrons we would see a much tighter relation. Such cleaning, 
though, is rather delicate and could introduce biases by itself. The trend is clearly 
there though even for the 'raw' data. 



maximum four-velocity is {'~iVz)max — 80. 

The theoretical model does of course not cover all details of the experi- 
ment. For example, in general the electrons also have velocity components 
parallel to the magnetic field; instead of making one-dimensional harmonic 
oscillations in the plane perpendicular to the current channel, the electrons 
will describe ellipsoidal paths. Fig. 16.11 shows the path of two electrons in 
the vicinity of an ion channel. But, overall, the electrons behave as ex- 
pected from the model considerations. Consequently, high-speed electrons 
are tightly coupled to the ion channels, as clearly illustrated by Fig. 16.1b . 

Figure 16.51 shows that the electrons are power-law-distributed at high 
energies, with index p = 2.7. The electrons at the high gamma cutoff are 
found where the ion current peaks, as may be seen from Fig. 16.41 The 
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maximum ion current is limited by the size of our box; larger values would 
probably be found if the merging of current channels could be followed further 
downstream. The PDF is not isotropic in any frame of reference because 
of the high anisotropy of the Weibel-generated electromagnetic field. The 
power-law in the electron PDF is dominant for 10 < 7 < 30. Likewise, a 
power law dominates the ion current channel strength, Jioni for 100 < Jjon < 
1000 (inset). A relation between the power-law distributions of these two 
quantities on their respective intervals is provided with Fig. 16.41 we see that 
the average four- velocity is proportional [straight line fit) to a power of the 
local ion current density on their respective relevant intervals, 10 < 7 < 30 
and 100 < Jjon < 1000. Their kinship stems from the fact that acceleration 
is local. The value J/on has a power-law tail, and its potential drives the 
high energy distribution of the electrons according to eq. 16.11 thus forming 
a power-law-distributed electron PDF. 

Measuring the rate at which the in-streaming ions transfer momentum to 
the ion population initially at rest allows us to make a crude estimate of the 
length scales over which the two-stream instability in the current experiment 
would saturate owing to ion thermalization. A reasonable estimate appears 
to be approximately 10 times the length of the current computational box, 
or about 1500 ion skin depths. Assuming that the shock propagates in an 
interstellar environment with a plasma density of ~ 10^ m~^, we may calcu- 
late a typical ion skin depth. Comparing this value with the upstream ion 
skin depth from our experiments, we find that the computational box corre- 
sponds to a scale of the order of 10^ m, or equivalently that the collisionless 
shock transition region of the current experiment corresponds to about 10® 
m. For an ion with a Lorentz factor 7 = 15, this length corresponds roughly 
to 40 ion gyro radii in the average strength of the generated magnetic field. 
But we stress that the in-streaming ions do not really gyrate, since they 
mainly travel inside the ion current channels where the magnetic field, by 
symmetry, is close to zero. Also, the strong electromagnetic fields generated 
by the Weibel instability and the nonthermal electron acceleration, which 
is crucial from the interpretation of GRB afterglow observations, emphasize 
the shortcoming of MHD in the context of collisionless shocks. 

By rescaling the experiment to physical units we find that the electrons 
are accelerated to maximum energies in the neighbourhood of 5 GeV. Even 
further acceleration may occur as ion channels keep growing down stream, 
outside of our computational box. 

The scaling estimates above depend, among other things, on plasma den- 
sities, the bulk Lorentz factor, and the mass ratio (mj/mg). A parameter 
study is necessary to explore these dependencies, but this is beyond the 
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Figure 6.5: The normalized electron particle distribution function downstream of 
the shock. The dot-dashed line is a power law fit to the non-thermal high energy 
tail, while the dashed curve is a Lorentz-boosted thermal electron population. 
The histogram is made from the four velocities of electrons in a thin slice in the z- 
direction of the computational box. The grey vertical line at = 15 indicates the 
injection momentum. The inset shows a similar histogram for ion current density 
sampled in each grid point in the same slice. The bump in the inset is a statistical 
fluctuation caused by a single ion channel. 



scope of this thesis. We thus stress that the extrapolations performed here 
are speculative and that unresolved physics could influence the late stages of 
the instability in new and interesting ways. 

When the in-streaming ions are fully thermalized, they can no longer sup- 
port the magnetic field structures. Thus, one might speculate that the radiat- 
i ng reg ion of the GRB afterglow is very thin, as suggested by Rossi and ReesI 



(200i). Furthermore, traditional synchrotron radiation theory does not ap 



ply to an intermittent magnetic field generated by the two-stream instability, 
since the electron gyroradii often are larger than the scales of the magnetic 
field structures. We em phasize the importance of the theory of jitter radia- 
tion (|Medvedevl . l2nnnl ). 
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6.4 Conclusions 

We have proposed a new acceleration mechanism for electrons in collision- 
less shocks. The theoretical considerations were suggested by particle-in-cell 
computer experiments, which also allowed quantitative comparisons with the 
theoretical predictions. We have shown that the nonthermal acceleration of 
electrons is directly related to the ion current channels in the shock transition 
zone. The results are applicable to interactions between relativistic outflows 
and the interstellar medium. Such relativistic outflows occur in GRB after- 



glows and in jets from compact objects (jPender et al.Ll2004^ . The suggested 



acceleration scenario mig ht overcome some of the problems pointed out by 



Baring and Brabvl ()2004f ) regarding the apparent contradiction between stan- 
dard Permi acceleration and spectral observations of GRBs. 

The mechanism has important implications for the way we understand 
and interpret observations of collisionless shocks: 

• The acceleration mechanism is capable of creating a power-law electron 
distribution in a collisionless shocked region. In the computer exper- 
iment presented here, a bulk flow with P = 15 results in a power-law 
slope p = 2.7 for the electron PDP. Additional experiments are needed 
to disentangle which parameters determine the exact value of the slope. 

• The acceleration is local; electrons are accelerated to a power law in 
situ. Therefore, the observed radiation field may be tied directly to 
the local conditions of the plasma and could be a strong handle on the 
physical processes. 



Our r esults strengthen the point already made by iPrederiksen et al. 



(200J) that the fractions of the bulk kinetic energy that go into in the 
electrons and the magnetic field, eg and e^, respectively, are not free 
and independent parameters of collisionless shock theory. Most likely, 
they represent interconnected parts of the same process. 

In the case of a weak or no upstream magnetic field, the Weibel-like 
two-stream instability is able to provide the necessary electromagnetic 
fields. We have shown here that the collisionless shocked region is 
relatively thin, and we suggest that the nonthermal radiation observed 
from GRB afterglows and relativistic jets in general is emitted from 
such a relatively thin shell. 



It is clear that the non-thermal electron acceleration, the ion current 
filament at ion, the magnetic field amplification/generation, and hence the 
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strong non-thermal radiation from the shock, is beyond the reach of MHD 
to explain. Whether the relativistic MHD jump conditions become valid on 
any larger scale is not possible to decide from the simulations presented in 
this chapter. 



Chapter 7 

Radiation from Collisionless 
Shocks 



In this chapter I investigate the radiation emitted from colhsionless plasma 
shocks. I construct a novel tool that generates radiation spectra directly from 
PIC simulations. This is an important step in order to link simulations with 
observations. Before utilizing this powerful tool, I perform numerous tests 
to support its credibility. I then investigate the nature of 3D jitter radiation 
in various setups, including snapshots from the PIC code. 



7.1 Introduction 

We first derive an expression for the frequency spectrum of the radiation 
emitted from an accele r ated relativistic charged particle. We initially follow 
Rvbicki and Lightman ( 19791 ). 



From Eq. 12.71 we have the total energy W emitted per unit solid angle 
per unit time from an accelerated charge 



dP d^W \RE{t)f 



(7.1) 



dfl dVtdt HqC 
The total energy W radiated per unit solid angle per unit frequency is then 

dildu hqc 

where 2tt comes from Parseval's theorem for Fourier transforms and the extra 
factor 2 comes from a symmetry argument of E(ci;) (no negative frequencies) 
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( Rvbicki and Lightman . 19791 ). E{uj) is the Fourier transform of E(t), which 
we define as 

E{uj) = — / E(t)e*"*rft. (7.3) 
Combining Eq. EH Eq. Oand Eq. Owe find that 

i?E(t)e*'^*cit 

n X [(n - /3) X /3] 



dfldui 



47r2 



167r3 



/3-n) 



(it 



ret 



(7.4) 



We now expand this into the retarded time frame: We need an expression 
for dt. From Section [2.3. II we recall that the retarded time is given by 



t' = t- R{t')/c. 
Taking the derivative of this expression with respect to t gives 



^ - 1 _ 

dt c dt 



^ _ 1 dR{t') dt' 
c dt' 'dt' 



(7.5) 



(7.6) 



We may obtain an expression for R{t') by expanding derivative of the identity 
i?2 = r2 into 2RR = 2R ■ R 
Thus kit') 



-2R ■ V since R 



ro 



R 



-ro 



-V. 



— n ■ V and we have 

dt={l-n- (3)dt'. 
With Eq. 17.51 and 17.71 we may write Eq. 17.41 as 

d^W 



dVtduj 



167r^ 



n X [(n - /3) X /3] ,^(t,„n 
(1-/3 -11)2 



ro(t')/c)^^/ 



(7.7) 



(7.8) 



In the exponential function, we have used t hat R{t') ~ |r| — n • fq. This 
is valid when R{t') ^ ro{t') () Jackson . 19991 ). which is always the case for 
observations at cosmological distances. The |r| can be neglected since in the 
exponential function, this corresponds only to an overall phase factor. 

In theory, this expression gives us the energy radiated per unit solid angle 
per unit frequency interval. However, two obstacles remain: We must know 
ro, (3 and as functions of time, and we must perform the time integration 
in Eq. 17.81 Regarding the first obstacle, we need an analytical expression 
for the particle orbit as a function of time. This may be found for a par- 
ticle moving in a homogeneous electromagnetic field in the non-relativistic 
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X Time 

Figure 7.1: The path of a charged particle moving in a homogenous magnetic 
field {left panel). The particle radiates a time dependent electric field. An observer 
situated at great distance along the n-vector sees the retarded electric field from 
the gyrating particle {right panel). As a result of relativistic beaming, the field is 
seen as pulses peaking when the particle moves directly towards the observer. 



limit. However, even here we are left with an integration that cannot be per- 
formed analytically. To make further progress, we assume that the external 
electromagnetic field is purely magnetic. 

The derivation from here to obtaining an expression for the e mitte d ra- 
diation as a function of frequency is rather extensive ( Jackson . Il999t ). In 
appendix El we follow this derivation with the extension of including the full 
angular dependency and we find the following result 







duodVL 


127r3 






duidVt 


127r3 



TiOp sin 6 



Ki 

3 



x/a/cos^ 



cos I 



cos( 



{cosOp 



3^2 



(7.9) 
(7.10) 



where 6 is the angle between n and the orbital plane 6^ = 2(1 — f3cos6), 
X = uvlO^/ {3c) and the gyro-radius '^mv/{qB). For [3^1 and 0-^0, 
this expression converges toward the so l ution one normally find in text books 
( Jacksonl . IT999[ iRvbicki and Lisfhtma^ . Il979l ). 

We restrict the solution to cos^ > 0. This is justified in the relativistic 
limit since the radiation from even mildly relativistic particles is beamed into 
a narrow cone around the velocity vector. Equation 17. 91 and Eq. 17.101 describe 
the synchrotron radiation spectrum as observed in the orbital plane. At low 
frequencies it has a characteristic power-law slope of uj = 2/3. Above some 
critical frequency it drops of as cue" 



-kuj 
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Figure 17.11 shows the electrical field E(t) from a particle that gyrates in 
a homogenous magnetic field as measured by an observer situated at great 
distance in the direction along n. 



7.2 Obtaining spectra from PIC simulations 

In general, the approximations above will not be entirely accurate. In the PIC 
simulations (and in the "real world"), the morphology of the electromagnetic 
field is often complicated, with a non- vanishing electric field. In this case, the 
paths of the particles are neither circular nor necessarily periodic. To obtain 
spectra from the simulations we need to take another approach. By sampling 
positions, velocities and accelerations from the simulation over a suitable time 
interval we can numerically evaluate the Fourier integral. First, we have to 
choose in which time frame we perform this integration: "proper" time t or 
retarded time t'. The corresponding integrals are given in Eq. 17.41 and Eq. 
17.81 At a first glance, Eq. 17.41 appears the most appealing since here we can 
use the standard Fast Fourier Transform (FFT) method, which is numerically 
faster. However, the term (1 — /3 • n) is numerically ill-behaved, and Eq. 17.41 
has a higher power dependency on this term than Eq. 17.81 Furthermore, for 
a relativistic particle moving towards the observer, all features in the emitted 
radiation appear strongly contracted in the t-frame and are thus considerably 
harder to resolve numerically. Finally, in order to use an FFT, the numerical 
representation of u needs to have uniform spacing, which is an inconvenient 
limitation. . 

Before we proceed, several other numerical issues need to be considered. 



7.2.1 Resolution constraints 



Several issues should be kept in mind when implementing the time integration 
of the retarded signal. Constraints exist that are tied to the uncertainty 
relation 

Au;At>l. (7.11) 



As explained by Rvbicki and Lightman ( 1979| ) "this uncertainty relation is 
not necessarily quantum in nature, but is a property of any wave theory of 
light" . As a result of this, the maximum frequency that can be represented 
in the resulting spectrum is limited by the number of points that one samples 
in a given time interval. This maximum frequency is known as the Nyq uist 



frequency and equals half the sampling frequency (e.g. iPress et al.lll986[ ). It 



corresponds to a rapid oscillating signal that repeats once per two sample 
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points in the time domain. Above the Nyquist frequency, the solution be- 
comes aliased, noisy, and in general invalid (see Fig. 17.21 and Fig. 17. 3^ . Thus, 
we set an upper limit on the frequency band to 

where At is the time between each sample point. The low end of the fre- 
quency band is also limited since the lowest frequency that can be represented 
is that of a wave that repeats only once in the whole time interval. In fact, we 
must sample for long enough to detect not only low frequencies in the signal, 
but also small differences between frequencies (the resolution). Thus, the 
length of time for which we sample the signal (Tg) determines our ability to 
resolve differences between frequencies. This defines the frequency resolution 

Au = ^. (7.13) 




Figure 7.2: The observed power spectrum from a single charged particle, gyrating 
in a magnetic field {left panel). Case A: Reducing the sample rate distorts the 
frequencies above the Nyquist frequency, marked by a vertical dotted line. Case 
B: If the total time interval becomes too small (smaller than the orbit period) the 
main frequency lobes become smeared out. The units on the axis are arbitrary. 



The effect of a too short total time interval can be seen from Fig. 17.21 Here 
we see that the individual spikes in the synchrotron spectrum are broadened 
out. So as a rule of thumb, the total length of the time sample should be at 
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least one over the lowest gyro-frequency in the particle ensemble for which 
we wish to obtain a spectrum. For an ensemble of, say, thermal particles this 
poses a lesser problem, since the individual frequency peaks become averaged 
out to a continuum anyway. 

10-'^ . . . H 



I 10-' 



10-" !,, ,, , , 

0.1 1.0 10.0 

Frequency 

Figure 7.3: An example of a numerical artifact introduced when trying to extend 
the spectrum to above the Nyquist frequency (dotted line). The effect is know 
as aliasing. The example spectrum is from electrons participating in the Weibel 
two-stream instability. 



7.2.2 Time domain windowing 

To determine the radiation power spectrum from an accelerated charged 
particle we sample the position, velocity and acceleration over a given time 
interval Ts = t2 — ti. However, even for periodic orbits we cannot presume 
that this time interval is an integer number times the orbital period. And in 
the general case of sampling many particles over the time interval, in an in- 
homogeneous electromagnetic field, there will most certainly be a high level 
of a non-periodicity. However, the theory behind Fourier transformations 
assumes that the signal being processed is a periodic waveform. Connect- 
ing the ends of non-periodic signal introduces a discontinuity in the signal. 
Such a discontinuity in the time-domain can only be represented by a super- 
position of all frequencies in the frequency domain. The result is that any 
peak frequency [main lobe) is smeared out over a frequency range and side 
lobes occur aroun d the main lobe. This is known as spectral leakage (e.g. 
Press et al. 198(t[ i. From experiments we have found that this leakage can 



alter the characteristic cu^^^ slope in the spectrum from a single synchrotron 
particle. 
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Figure 7.4: Windowing of a signal sampled over a time block ensures that the 
signal becomes periodic and continuous in the end points. This greatly reduces 
the spectral leakage. 



To solve this problem, we apply a filtering function on the time dependent 
signal, to ensure that the signal becomes exactly zero in the end points of 
the chosen time sample block (and thus periodic), so that they match up and 
the transition is continuous (Fig. I7.4|) . The technique is known as windowing 
and is a standard tool in discrete signal analysis. The choice of window 
function depends on the problem. Several standard functions exist (Hanning, 
Hamming, Kaiser-Bessel...). After trying several of these, as well as other 
functions on the synchrotron radiation case we have found a good ad hoc 
window function 



W{t) = Exp 



u-u 



(7.14) 



with the parameters /c = 3 and m = 6. 

After applying a window function to the signal it is important to renor- 
malize the sample in order to retain the amplitude of the spectrum. Thus 
we multiply the final spectrum with the factor where 

Aw = 4^^. (7.15) 

Cw{t)dt 



With the window function applied on the sampled data, the expected 
coi2/3 synchrotron slope is effectively restored. 
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Figure 7.5: To obtain a spectrum from the PIC simulation at time t, particles are 
traced over a time interval around t with high temporal resolution. Positions, ve- 
locities, and accelerations are sampled and from eauation l7.8l we find the radiation 
spectrum. 



7.2.3 Implementation and testing 

The procedure above has been implemented in two steps. The first part is 
to sample position, velocity and acceleration for each particle. In general, 
one could dump particle data directly from the simulations. This, however, 
requires not only vast amount of disk storage, but one also needs to run 
the simulations with a time step much smaller than otherwise needed. The 
reason for this is discussed in Section 17.2.11 Typically, the time step in the 
spectrum tracer is of the order of 100 times smaller than in the PIC code. 
Since the simulations already take several weeks to run, we choose to take 
another approach (Fig. 17. 5j) : When a simulation has finished, we choose 
a snapshot for which we want to obtain the radiation spectrum. Loading 
fields and particles from the chosen time step, we trace a given number of 
the particles back and forth in a time segment with very high temporal 
resolution, while keeping the electromagnetic fields frozen. In this sense, 
the particles can be seen as test particles. At each small time step we then 
store the particle positions, velocities and accelerations in a file. This data is 
read by another program that evaluates the Fourier integral (in the retarded 
time frame) given a direction n towards the observer. Finally we add the 
spectrum of each particle into a total spectrum. Adding the spectra from 
each particle linearly is valid as long as the phase of each contribution is 
completely uncorrelated to the others. This is clearly an assumption, which 
would be violated for example in a maser arrangement. The spectrum may 
finally be Lorentz/Dobbler boosted to any observer rest frame. 
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Figure 7.6: The spectrum obtained from a single particle gyrating in a magnetic 
field from the code (blue). The thick red line shows the theoretical synchrotron 
spectrum. The dashed line indicates the power-law slope 2/3. The inset in the 
upper left corner shows that the spectrum consists of a discrete set of spikes that 
are integer overtones of the gyrofrequency, in this case ujb = 0.094 (arbitrary units) 
{dotted vertical lines). The particle has a Lorentz factor 7 = 5. All numbers are 
in simulations units. 



Test: Synchrotron radiation 

To verify the ability of radiation code to produce correct spectra from ac- 
celerated charges, a series of tests have been performed. The first test is 
that of synchrotron radiation. We place a single relativistic particle in a 
homogeneous magnetic field (Fig. 17. The theory predicts that the emit- 
ted energy spectrum should consist of a collection of peaks, positioned at 
the gyro-frequency and higher harmonics. The amplitude of the individual 
peaks follow the continuum profile given by Eq. 17.101 Fig. 17.61 shows the 
spectrum we obtain from the code for a single relativistic charged particle 
in a homogenous magnetic field. We find very good agreement between the 
spectrum obtained from the simulation (blue line) and the theoretical syn- 
chrotron spectrum from Eq. 17.101 (red line). As expected, the peaks are found 
at harmonics of the gyro eigenfrequency and the low energy part follows a 
power-law slope u'^^^. 

If one observes not a single synchrotron particle, but rather a whole en- 
semble of particles, the characteristics of the spectrum change. If all the 
particles in the ensemble have the same energy, but have momentum vectors 
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that are isotropic distributed, the characteristic low-energy power-law slope 
of the integrated spectrum is altered. Integrating Eq. 17.91 and 17.101 over all 
pitch angles yields a spectrum with a low-en ergy part of the spectrum tha t 
has a power-law slope of 1/3 rather than 2/3 (jRvbicki and Lightmanl . ir979| ). 
The results from the simulations are in excellent agreements with this. Figure 
17.71 shows the spectrum from 512 particles, sampled from a mono-energetic 
distribution function (7 = 10). 




0.01 0.10 1.00 10.00 

Frequency 



Figure 7.7: The low energy part of the spectrum from an ensemble of electrons 
with 7 = 10 and random orientation of the momentum vector. The ensemble is 
placed in a homogeneous magnetic field. The dotted line indicates a power law 
slope of 1/3. 512 particles were traced for this spectrum. 



Test: Bremsstrahlung 

Another important test that our radiation tool must pass is to reproduce 
the theoretically predicted spectrum of bremsstrahlung. Bremsstrahlung is 
emitted when the path of a charged particle is perturbed by an electric field 
(e.g. an ion). In the present version of the PIC code (and in most other PIC 
codes), the particles cannot interact directly but only through the grid. For 
simulations where each cell in the computational domain is quasi-neutral, 
bremsstrahlung is not included. The next generation of out PIC code will 
be able to include particle-particle interactions to some extent (see Chapter 
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IHl). We can, however, still test the radiation module by placing a single 
ion/electron pair in a simulation box. 

Where synchrotron radiation is a periodic phenomenon tied to a magnetic 
field, bremsstrahlung is non-periodic and connected to electric defiections. 
So in a sense, this test complements the synchrotron test very nicely. Figure 




Time 



Figure 7.8: A relativistic electron passes an ion (left frame). The electrical far- 
field peaks as the electron is accelerated by the ion (right frame). 



17.81 shows the example setup where an electron passes an ion and becomes 
subject to acceleration. In the relativistic limit the theoretical spectrum for 
this process may be derived using the beautiful theory of virtual photons 
( Rvbicki and Lightman . 1979[ ): In the rest frame of the electron, the electric 
field from the ion appears as a short, Lorentz contracted, electromagnetic 
pulse that can be view as a virtual photon. This virtual photon Compton 
scatters of the electron to produce the emitted radiation. Transforming back 
to the lab frame of the ion gives the relativistic bremsstrahlung emission from 
the electron 



(7.16) 



Here, Z is the number of charges in the ion, h is the minimum distance 
between that ion and the electron and Ki is the modified Bessel function 
of the second kind. We note that a factor of Ki is missing in Eq. 5.23 
on page 164 in Rvbicki and Lightman ( 1979| ). In Eq. 17.161 the _L indicates 
that we have neglected the contribution from the component of the ions 
electric field that is parallel to the path of the incoming electron. This 
approximation is valid in the ultra relativistic limit (7 ^ 1). Figure 17.91 
shows the simulation results {thick blue) compared with Eq. 17.161 {thin red) 
for different electron energies (7^=0.5, 2.5, 5.0, 7.5, 10.0). We find good 
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agreement between simulations and the theoretical result for > 1. The 
deviation for 7 — > 1 is because of the approximation in Eq. 17. 161 above. 
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Figure 7.9: The theoretical {thin red) and simulated (thick blue) bremsstrahlung 
spectrum from an electron passing an ion. Four cases have been tested, differing 
only in initial electron energy (relativistic momentum ^v=0.5, 2.5, 5, 7.5, 10). 



Test: Undulator radiation 



A final but crucial test is whether the radiation code can reproduce the 
spectrum from undulator radiation. Undulator radiation is emitted when 
a relativistic electron traverses a periodic magnetic field so weak that the 
deflection angle stays within the relativistic beaming cone. The phenomenon 
is well known, and widely used, in synchrotron storage ring experiments. 
Figure 1 7. 101 shows a sc hematic view of the generation of undulator radiation 



( Attwood etin I1993I ). A full introduction to the field of undulator and 



wiggler radiation is given bv lAttwoodI (l2nnn> ). 

The spectrum from a relativistic electron moving in a transverse periodic 
undulator field (with magnetic wavelength A„ and maximum strength Bq), is 
characteri zed by a well c onfined peak near Uu, determined by the undulator 
equation ( Kincaidl . Il977 ) 



2c7' 



A„(l + i^72 + 72 



(7.17) 
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Figure 7.10: A schematic view of a undulator experiment in a cylotron. A cluster 
of magnets are aligned along the path of the electron. When the electron passes 



the ar ray of magnetic, it radiates at the undulator frequency. From lAttwood et al 



Here, 6^ is the angle between the path of the electron and the line of sight 
and K is the periodic magnet parameter, defined as 

K = XuqBo/{27rmec). (7.18) 

In the weak-field case [K < 1) the deflection angle of the undulation is 
smaller than the beaming opening angle I/7 and the radiation is confined 
into the angle I/7. For K > 1 the radiation (called wiggler radiation) is 
emitted into a cone of half-angle K/'-f. In this section we focus on K < 1. 
The width of the peak is determined by how many periods TV the electron 
path is sampled over. The more undulations, the better defined the peak 
becomes. 

Figure 17. 1 II shows the result of simulations where a single, relativistic elec- 
tron (7 = 10) moves in a periodic magnetic field. The undulator wavelength 
is Xu = 11 grid-zones and the periodic magnet parameter K = 0.0015 (Eq. 
I7.18p . The four panels in the figure represent different angles between the 
electron path an d the line of sight (6'„=0, 0.01, 0.05, and 0.1). As expected 
( AttwoodlboOOi l the spectrum peaks at the frequency given by Eq. l7.17l (dot- 
ted lines). As 6u is increased, higher harmonics of Uu becomes visible. We 
thus find excellent agreement between the undulator equation Eq. 17.171 and 
the simulations. 
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Figure 7.11: Undulator spectra generated when a 7 = 10 electron traverses a 
periodic magnetic field. The spectra are obtained from simulations with different 
viewing angles 9u= 0, 0.01, 0.05, and 0.1 (top to bottom panel). The dotted lines 
show the expected undulator frequency lOu and its higher harmonics. 



7.3 Jitter radiation 

After having tested the radiation generation module thoroughly, we now focus 
on some more scientific issues. The PIC plasma simulation code, in combina- 
tion with the radiation generation module, provides us with a powerful tool 
for testing various non-linear problems. 

We choose to address several problems, all connected to the theory of 3D 
jitter-radiation from GRB plasma shocks. Ji tter radiation was introduced 
to the GRB community by Medvedev ( 200Cl( ). It covers the regime where 



the magnetic field is inhomogeneous on scales smaller than the Larmor ra- 
dius (see Fig. I7.12j) and the electron's transverse deflections in these fields 
are much smaller than the relativistic beaming angle. The foundation and 
motivation for such a theory was based on the small-scale nature of the mag- 
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netic fields generate d by the two-stre am instability. In the one-dimensional 
analytical approach, Medvedev ( 200Cl( ) found that the low frequency slope of 
the spectrum could be steeper than the 1/3 slope of synchrotron radiation. 
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Figure 7.12: An example of the complicated paths of the electrons under the 
influence of electromagnetic field generated by the Weibel two-stream instability. 



In this section, we first perform a parameter study to investigate the 
nature of the radiation emitted from an ensemble of relativistic electrons in 
different small-scale magnetic configurations (all with K < 1). 

Secondly, we compare the radiation that is emitted from a electron pop- 
ulation in a generic, turbulent magnetic field (power-law distributed power- 
spectrum) with a two-stream generated magnetic field obtained from a snap- 
shot from the PIC simulations. The question is: Is it reasonable to apply 
the theory of jitter-radiation directly on GRB shocks, assuming some generic 
structure of the magnetic field or must one include the complicated spatial 
details of a two-stream generated magnetic field? Our concern is that the 
ordered sub-structures in the two-stream generated magnetic field may affect 
the radiation spectrum. 



7.3.1 The 3D jitter spectrum 

Jitter radiation is in family with undulator radiation, except that the mag- 
netic field has power in many Fourier harmonics rather than only one, and in 
general follows a power-law in Fourier space. We define the power-spectrum 
of the magnetic field 

PBik) = Csk^ (7.19) 

(note that iMedvedevI (j2nn(l l defines the spectrum as Bk = Csk^, which 
means that there is a factor 2 difference on fi compared to our definition). 
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Figure 7.13: 2D-slices of a turbulent magnetic field for three different values of 
the spectral power law index /i: /i = —2 (red noise, top), /U = (white noise, 
middle) and /i = 2 (blue noise, bottom). 



In the original calculations. iMedvedev (l2000f) assumed th at /i > 1. However, 
large-scale PlC-s„„ulations7IiniSStoS^ ») have shown that 
the magnetic field generated by the non-linear Weibel two-stream instability 
is indeed power-law distributed but has /i < 0. 

In this section we investigate the radiation from an ensemble of relativistic 
electrons, moving in a random magnetic field, for different values of the 
spectral power law index fi (see Fig. 17.1^1) . The result is a three-dimensional 
generalization of the jitter-theory that differs substantially from the one- 
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Figure 7.14: The jitter spectra from a mono-energetic ensemble of electrons (7 = 
3) moving in a turbulent magnetic field. The different graphs are for simulations 
with different values of /i. The amplitude of the magnetic field is set so that the 
periodic magnet parameter K < 1 for all /c-nodes (Eg. I7.1(j|) . 



dimensional result found by Medvedev (200y). In this thesis, we do not 
i nclude the regime where an additional large-scale amplitude field is present 
( Medvedevl .l2nnf]). 

We have performed simulations of an isotropic mono-energetic ensemble of 
electrons (7 = 5). The isotropic distributions may be seen as an integration 
over all angles of the undulator equation. The electrons are placed in a 
simulation box with 200 x 200 x 800 grid-zones with a turbulent magnetic field 
that has a Fourier distribution as in Eq. 17.191 Seven different values of /i have 
been examined. The amplitude of the magnetic field is set so that the mean 
magnetic field energy is the same in each simulation. The periodic magnet 
parameter K is less than 1 for all Fourier nodes. The range of the power-law 
is limited to kmin < k < k^ax where kmin and k^ax correspond to minimum 
and maximum frequencies that can be represented by the simulation grid 
(Fig.mni). 

The resulting 3D jitter radiation spectrum from the electrons can be seen 
in Fig. 17.141 Quite interestingly, all the spectra have similarity with the 
bremsstrahlung spectrum for low energies (fiat spectrum). At higher ener- 
gies, the fi > cases drop off very sharply - faster than the bremsstrahlung 
spectrum. For fi < the high-energy part of the spectrum is found empir- 
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Figure 7.15: Increasing the magnetic field increases the magnet parameter K 
(Eq. I7.18jl . The left panel shows how the amplitude of the undulator spectrum 
increases with {B'^). When K becomes larger than 1, a transition to the wiggler 
radiation spectrum occurs. The right panel shows the same spectra divided with 



ically to be a power-law with a slope of uj^^^. The reason why the /i > 
peaks at so much higher frequencies than for yU < can be understood from 
Eq. 17.171 which shows that the characteristic frequency scales as k (~ l/A.^). 

We have also tested how the spectrum depends on the magnetic field am- 
plitude. From Fig. l7.15l we find that the amplitude scales with as expected 
from Larmor's radiation formula. Above the limit where the deflection angle 
becomes comparable to, or larger than the beaming angle {K > 1), the spec- 
trum enters the wiggler- do main. The wigg ler radiation spectrum is close to 
the synchrotron spectrum ( Attwoodl . 2000[ ). 



We finally test the 7-dependency. Figure l7.16l and Fig. l7.l7l show the jitter 
radiation spectra for /i < and /i > for electron ensembles with different 
energies. Even though the spectra in the two cases differ, they both scale the 
same way: The spectra are shifted with 7^ in frequency in agreement with 
Eq.lOKand Ea. lTTTjl . 



7.3.2 Jitter radiation from collisionless plasma shocks? 

In this section, we compare the radiation spectrum from an ensemble of elec- 
trons, moving in two different magnetic topologies: 1) the magnetic field 
genera ted by the Weibel two-stre am r nagnetic field Bn (see S ection El in this 
thesis, Medvedev and Loebl 19991 and Frederiksen et al. hoO^ and 2) a ran- 
domized magnetic field Bi that has the same average energy density and 
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Figure 7.16: The radiation spectrum from fom' different mono-energetic ensem- 
bles of electrons (7=8, 4, 5, and 9) with isotropic distributed momentum-vectors, 
placed in a turbulent magnetic field with spectral slope fj, = —3 {left panel). The 
right panel shows the same four graphs but shifted in frequency with a factor 7"^ 
(renormalized) . 
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Figure 7.17: The radiation spectrum from four different mono-energetic ensem- 
bles of electrons (7=8, 4, 5, and 7) with isotropic distributed momentum vectors, 
placed in a turbulent magnetic field with spectral slope fj, = +3 {left panel). The 
right panel shows the same four graphs but shifted in frequency with a factor 7"^^ 
(renormalized). 



spectral power spectrum as the two-stream generated field. The purpose of 
this exercise is to test if there are any differences in the resulting radiation 
spectra. Bq is readily available from the PIC simulations. Constructing Bi 
is more tricky. To make sure that Bi has exactly the same average energy 
density and power spectrum properties as the PIC-field Bq, we generate the 



94 



CHAPTER 7. RADIATION FROM COLLISIONLESS SHOCKS 



random field in the following way: 

The basic idea is to Fourier transform the magnetic field from Bq, make a 
random phase shift of all fc-nodes in the Fourier domain and then transform 
back to real-space. This is easily done but has one serious flaw: we cannot 
be sure that the resulting magnetic field obeys V ■ Bi = even though Bq 
does. The solution is to work on the vector potential rather than directly on 
Bq. The vector potential Aq is a vector field that satisfies Bq = V x Aq. 
Inverting this equation requires Bq to be periodic. In the PIC simulations, 
Bo is already periodic in the two-dimensions perpendicular to the jet flow 
{x,y). To make the magnetic field periodic in the 2;-direction we apply a 
windowing function that connects the magnetic field-lines through the z- 
boundary. Away from the boundary, the magnetic field remains unchanged. 
We have tested that this filter does not alter the power spectrum of Bq except 
for a few percent in the very highest frequencies. Ensured that Bq is periodic, 
the Fourier transform of Aq into Aq is 

Ao(k) = FFT(Ao) = ik x FFT(Bo)/|kp. (7.20) 

k is the wave vector and FFT represent a fast Fourier transformation. 

We then scramble the vector potential by a random phase shift of all the 
Fourier harmonics 

Ai(k) = Ao(k)e*2'^k/ (721) 

where / is a matrix of random numbers in the range [0; 1]. The size of / 
matches the total number of data points in the three-dimensional Bg-array. 
Bi is then given as Bi = V x Ai. The transformation is performed in IDL 
(Interactive Data Language by Research Systems Inc.) and this adds a few 
complications to the process, mainly connected to the representation and 
renormalization of the Fourier transforms. To ensure that no errors have 
been introduced in the process of obtaining Bi, we check and find that 

• Bq and Bi have identical Fourier transforms and equivalently, identical 
spectral power distribution, 

• Bi is real if Bq is real, 

• the magnetic energy is conserved through the transformation J ^IdV = 
J BldV, and 

• Bi is divergence free (V ■ Bi = 0). 

Figure 17. 1 81 shows the input magnetic field Bq and the resulting Bi. Three 
slices are shown for each field, sampled at different depths in the shock. 
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Magnetic energy 




Figure 7.18: The spatial distribution of magnetic energy in three sUces sampled 
at different depth in the plasma shock. The top panel shows Bq generated by the 
two-stream instability. The lower panel shows a random magnetic field Bi with 
the same spectral power distribution as Bq. 



Figure l7.19l shows the power spectrum of Bq and Bi. Clearly, the randomized 
magnetic field Bi has the same statistical characteristics as the magnetic field 
from the simulations. 

In Fig. I7.2UI we compare the radiation spectrum from a thermal electron 
ensemble in Bq, with the spectrum of the same ensemble in Bi. Bq is from 
PIC simulations of a plasma shock (F = 3 simulation from Chapter E)). Bi 
is created as described above. Both spectra are generated by tracing 20000 
electrons. To make the spectra as realistic as possible, the electron positions 
and velocities are sampled from a snapshot of the PIC simulation. 

We find that the two spectra peak at the same frequency but diverge 
considerably for high frequencies where the radiation spectrum from the ran- 
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Figure 7.19: Power spectrum of the magnetic field generated by the Weibel two- 
stream instability Bq {full line) and the random phase shifted field Bi {dotted 
line). 




Figure 7.20: Radiation from PIC simulations {solid blue) compared to radiation 
from a random field with the same statistical properties {dashed red). Both spectra 
are generated by using the same ensemble of 20000 electrons. 



dom field Bi has a harder cutoff. The reason for this is because of the higher 
peak- values in the magnetic field from the PIC simulations (see Fig. I7.18|) . 
Even though Bq and Bi have the exact same power-spectrum, there exists 
phase-correlations in Bq that give the field a more coherent structure (Bi 
has a larger volume filling factor). We know also that the frequency where 
the synchrotron spectrum peaks Uc scales linearly with the magnetic field 
strength. Finally, the high-energy electrons are primarily found near ion- 
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current channels where the magnetic field peaks (see Fig. 13.11 16.11 and 16. 4j) . 
We conclude that it would appear that one may use a random magnetic field 
for analytical calculations, but only as a first approximation. We also stress 
that the simulation results we have used are from simulations that do not 
include the full shock ramp. The spectra may diverge even more for radiation 
from a full shock. 

Radiation due to Electric Fields 

Electric fields are most often neglected in generation of radiation from as- 
trophysical shocks. In PIC simulations we have found that for coUisionless 
shocks where the magnetic field is generated by the Weibel two-stream insta- 
bility, electric fields exist with an energy density of the order of 10% of the 
generated magnetic field. Figure 17.211 shows two radiation spectra emitted 
from a coUisionless shock. The electron distribution and electromagnetic field 
that lie behind the radiation are taken directly from self-consistent PIC code 
simulations. The left panel of 17. 2 II shows a spectrum where only the magnetic 
field is taken into account, while the right panel shows the spectrum when 
the electric field is also included. For a pure magnetic field we immediately 
identify the wiggler/synchrotron signature from the oj^/'^ low energy slope. 
The main effect from including the electric field is that the low energy slope 
is flattened to ~ uj^^^ and the high-energy exponential cut-off is softer and 
goes to higher energies. From these plots we conclude that it is important to 
also include the electric fleld generated in the Weibel two-stream instability. 

7.4 Spectra from the PIC simulations rescaled 
into 'real space' 

The ultimate goal is to obtain spectra from the PIC shock simulations that 
can be compared directly to observations. Even though it is beyond the scope 
of this thesis, this will eventually allow us to put constraints on the physics 
of GRB afterglows. 

In the PIC simulations, several parameters have been scaled to simulation 
units. To make the simulation results more directly comparable to real ob- 
servations, it is necessary to correct for this. Rescaling the radiation spectra 
obtained from simulations involves several effects: 

Simulation time-scales: To rescale the time-scale of the simulations to 
real space, we must divide all frequencies with the plasma frequency 
in the upstream part of the simulation box uj' and multiply with the 
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Frequency Frequency 

Figure 7.21: The spectrum of a mono-energetic (7 = 3), isotropic electron popu- 
lation in the electromagnetic field generated by the Weibel two-stream instability. 
In the left panel, only the magnetic field is included whereas the right panel in- 
cludes both the magnetic and electric field. The grey shaded line in the left panel 
is a copy of the spectrum in the right panel and visa verse. 20000 particles were 
traced in each spectrum. Both axis are in arbitrary units. 



real space (ISM) plasma frequency ujp. That is, we have to multi- 
ply with a factor Ri = Up/u'p. In the non-linear phase of the Weibel 
two-stream instability, the ions are dominating the evolution of the 
instability. Therefore, we set u'^ equal to the upstream ion plasma fre- 
quency from the simulations {u'p = 0.025). Since we are focusing on 
the afterglow domain in the simulations, we choose the ion plasma fre- 
quency of the inter-stellar medium for the real space plasma frequency, 
{ujp - 1300 s-^). In this case, Ri = 1300/0.025 ~ 5 ■ 10^ 

Ion to electron mass ratio: For numerical convenience, the ion (proton) 
to electron mass ratio are rescaled (see Section |2I). The ion/electron 
charge ratio remains unity. From the simulations it is found that 
the ions dominate the Weibel two-stream instability. Therefore, we 
can see the mass rescaling as an unphysical rescaling of the electron 
mass/charge ratio. In the code, the ion/electron mass ratio is typically 
16. Using the real mass ration 1836 would mean that the electrons 
become accelerated to higher Lorentz factors (but the same energy) by 
a factor 1836/16. In Section [7.3.1l we found that the peak frequency in 
the spectra scales with 7^. Thi s is consistent with synchrotron radiation 
( Rvbicki and Lightman . 1979[ l: In synchrotron radiation, the peak fre- 



quency is proportional to 'j'^/rrig. Thus we must effectively shift the ra- 
diation spectrum up in frequency by a factor R2 = (1836/16)^ = 1.5-10^ 
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Relativistic Doppler shift: All the simulations are carried out in the rest 
frame of the shock, moving with a bulk lorentz factor F towards the ob- 
server (only GRB jets moving directly towards us will trigger the space 
telescopes as a result of the relativistic beaming). We must therefore 
shift the simulation radiation spectra with a relativistic Doppler cor- 
rection term ^ 2r. 

In summary, to compare the spectrum we have obtained from the simu- 
lations of a r = 15 shock propagating through the ISM, we must shift the 
frequency axis with a factor R = R1R2R3 — 2-10^^. This puts us somewhere 
near the optical frequency band. In Fig. 17.221 we show the radiated spec- 
trum, sampled from 20.000 electrons. After applying the rescaling, we find 
that the spectrum peaks in the far infrared area. Below the peak, we find 
a power-law segment with slope P{uj) oc cj^/^. This is interesting because it 
is steeper than the synchrotron 1/3 slope and thus might shed light on the 
issue about "the line of death" ()Preece et al. . 1998| ) for many bursts. For 



frequencies above the peak frequency, the spectrum continues into the near 
infrared/optical band, following a power-law with slope P{uj) oc uj~^ with 
/3 - 0.7. 

From the power-law slope in the spectrum with {j3 = 0.7), the standard 
procedure in the GRB community would find that the electrons have a power- 
law distribution, with slope p = 2.4 (determined by solving [p — l)/2 = (3). 
The standard conclusion would be that this is consistent with Fermi acceler- 
ation. However, we strongly emphasize that the electrons that produced the 
spectrum in Fig. 17.221 are not Fermi accelerated. The electrons are instan- 
taneously accelerated and decelerated in the highly complicated electric and 
magnetic field near the ion current channels under emission of strong radi- 
ation (see an example of the complicated electron-paths in Fig. I7.12|) . This 
differs substantially from the iterative acceleration in Fermi acceleration. 

We note that the box of the PIC simulation behind the spectrum pre- 
sented in Fig. 17.221 is not large enough to cover the whole shock ramp (ions 
are not fully thermalized when they leave the simulation domain). This 
means that additional features might be added to the spectrum for even 
larger simulations. Moreover, the spectrum does not include synchrotron 
self-absorbtion. 



7.5 Summary 

I have developed, and tested, a tool can create radiation spectra directly 
from particle-in-cell simulations. By tracing particles, the code performs the 
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Figure 7.22: The spectrum from a F = 15 plasma shock propagating into an 
ISM-like medium. The frequency scale is boosted in accordance with the physical 
arguments given in the text. The spectrum peaks in the far infrared. Below the 
peak, we find a power-law segment with slope P{uj) oc w^/^. For frequencies above 
the peak frequency, the spectrum continues into the near infrared/optical band, 
following a power-law with slope P{ijj) oc w"^, where f3 = 0.7. 



exact radiation field Fourier integration (Eq. 17. 8j) without having to make 
assumptions about the magnetic field, particle orbit, beaming angle, and so 
forth. 

The radiation tool has been thoroughly tested and successfully reproduces 
synchrotron radiation, bremsstrahlung, and undulator radiation from small- 
angle deflections. 

I have used the tool to investigate the properties of three-dimensional 
jitter radiation in magnetic fields with different turbulent configurations. 
By tracing isotropic momentum distributions of electrons through a random 
magnetic field with a power spectrum that follows a power-law distribution 
in the Fourier domain Pb(/c) oc k^, 1 have determined the resulting spectrum. 
I have focused on the weak field limit (where the deflection angle is less than 
the beaming cone angle, K < 1). For all values of /i, the radiation spec- 
trum has a flat low-energy part (similar to bremsstrahlung). For /i < 0, the 
high-energy part of the spectrum follows a power-law with a slope a ~ /i — 1 
independent of the electron energy. For /x > 0, the flat spectrum continues to 
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higher frequencies than for jj < 0. The spectrum has a hard cut-off for high 
frequencies. With the present simulations, it is not possible to determine if 
the spectrum has peaks for even higher frequencies. 

For all values of fi and K < 1, the spectrum shifts in frequency with 
7^ and scales with the amplitude of the magnetic field squared B^, in good 
agreement with Larmor's formula for radiated power. In the limit of large 
deflections [K > 1) the spectrum eventually converges through the wiggler 
spectrum to the ordinary synchrotron spectrum. 

For all the spectra presented in this chapter, I have checked for conver- 
gence: For each spectrum obtained I have made a second run with twice 
the number of particles. If the two spectra have different shape, the pro- 
cess is repeated by doubling the particle number until the result converges. 
This process is important since radiation from high-energy particles is highly 
beamed. This is prone to distort the high-energy radiation spectrum because 
of bad statistic. 

Computationally, the generation of spectra from simulations is not an 
easy task. One may say that it is an example where nature is not nice 
to scientists. Especially, obtaining a radiation spectrum from high-energy 
particles is not easy since it is necessary to integrate for a time period that 
grows linearly with 7 in order to sample at least a couple of orbits. But at 
the same time it is crucial to sample a very high number of time-steps per 
time unit, to resolve the peak in the electrical field, which is very narrow 
(7"^) because of relativistic beaming. So in total, the computational cost 
grows with 7^^. The good news is that each particle and frequency bin may 
be treated independently, which makes the problem embarrassingly parallel 
(in the jargon of High Performance Computing). 

The method of tracing particles in an electromagnetic field that is fixed 
to a single snapshot is clearly an approximation. For high frequency spectra, 
this is not a problem, since we sample for a very short period of time. For 
low frequencies, it is a potentially more severe approximation since it doesn't 
make sense to trace the particle for too long time in a field that would 
have changed significantly during this process. Improvements can be made 
in several ways. One method is to pick out a number of particles before 
the large simulation starts and then integrate these particles with a high 
number of sub-steps. Another approach would be to interpolate field values 
to the particles position not only in the spatial dimensions (which is already 
implemented) but also in time from different field snap shots. 

An interesting step for the future is to determine the polarization in the 
synthetic spectra from the PIC simulations. Even though time did not allow 
for it in this thesis, it should not be difficult with the tools that are already 
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developed. One could then investigate how the radiation spectrum, including 
polarization change as functions of the viewing angle measured from the jet 
propagation direction. 

In the near future, the computational resources will reach a level where 
one can resolve fully three-dimensional collisionless shocks. It then becomes 
possible to make composite spectra from different simulations with varying 
bulk Lorentz factor as a function of viewing angle. This will allow us to 
test many interesting aspects with regard to jet structure and polarization 
predictions from shock-generated electromagnetic fields. 

I end this chapter by emphasizing that strong magnetic field generation, 
particle acceleration and emission of non-thermal radiation are unavoidable 
consequences of the collision of two, initially unmagnetized, plasma shells. 



Chapter 8 

A next generation PIC code 



8.1 Introduction 

Over the last couple of years the Copenhagen group has been using PIC mod- 
els that include electromagnetic fields and ch arged particles to understand th e 
plasma microphysics of coUisionless shocks Frederiksen et al. ( 20021 2004) 



plasma micropnysics qi comsioniess snocKs i.f reaeriKsen et ai.l ijzuu/i I/UU4J I: 
Hededal et all (|2004l ): iHededal and Nishikawal (|200,^ . It has turned out to 



be a very successful tool, but it is still limited in the scope of phenomena 
that can be addressed. Even though a large class of astrophysical environ- 
ments are indeed coUisionless, scattering and collision processes do play an 
important role in several key scenarios. Examples are given below. Another 
key ingredient, which has been missing in charged particle simulations, is 
a full treatment of photon propagation. It can be argued that photons are 
represented directly on the mesh by electromagnetic waves, which certainly 
is correct. But the mesh can only represent waves with frequencies smaller 
than the Nyquist frequency. The physical length of a typical cell has in our 
applications typically been 10^ — 10^ cm and hence it is clear that only low 
frequency radio waves can be represented. High-frequency photons have to 
be implemented as particles that propagate through the box and interact, 
either indirectly through messenger fields on the mesh, or directly with other 
particles. A valuable consequence of modelling the detailed photon trans- 
port is that extraction of electromagnetic spectra is trivial. Even in cases 
where the photon field is only a passive participant, this fact should not be 

underestimated as it enables direct comparison with obse rvations. 

There exists Monte Carlo based particle codes (see e.g. IStern et all t99± 



and references therein) that address various particle interactions, but one 
of their main shortcomings is the poor spatial resolution. This makes it 
impossible to couple the particle aspects to a self-consistent evolution of the 
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plasma. 

Our goal has been to develop a framework where both electromagnetic 
fields and scattering processes are included in a consistent way. We can 
then correctly model the plasma physics and the radiative dynamics. The 
scattering processes include, but are not limited to, simple particle-particle 
scattering, decay and annihilation/creation processes. Our new code is not 
limited in any way to charged particles, but can also include neutrals such 
as photons and neutrons. 

In the next section we describe some of the physics that can be addressed 
with this new code. In section 18.21 we discuss how the code has been imple- 
mented, the general framework and in detail, which physical processes that 
are currently implemented. In section IH^ we present the results of a prelim- 
inary toy experiment that we have performed to validate the code. In the 
last section IH3 we summarize. 



8.1.1 Motivation 

Before we continue and describe in detail the methods, physics and test 
problems we have implemented and used, it is important to consider the 
general class of scenarios we have had in mind as motivation for developing 
the code. There are several key objects, where only the bulk dynamics is 
understood, and we are lacking detailed understanding of the microphysics. 



Internal shocks in Gamma-Ray Bursts 



In the internal/external GRB shock model, the burst of gamma-rays is be- 
lieved to be gene rated when relativis t ic she ll s collide and dissipa t e thei r rela- 
tive bulk energv lRees and MeszarosI ( 1992| ): Meszaros and Reei f 1993[ ). The 
nature of the radiation is presumably inverse Compton scattering and syn- 
chrotron radiation. Particle/photon interactions mi ght also play an impor- 
tant r ol e in the very early a fterglow as suggested bv [Thompson and Madau 
( 2000l ): lBeloborodo\] J 2OO2I ): Even though the medium that surrounds the 
burst (ISM or wind) is optically very thin to gamma-rays, a tiny fraction of 
the gamma-rays will Compton scatter on the surrounding plasma particles. 
This opens up for the possibility of pair-creation between back scattered 
and outgoing gamma-rays. The creation of p airs niay inc rease the rate of 
back scattered photons in a run-away process ISternI (l2003t ). The Compton 
scattering may accelerate the pair-plasma through the surrounding medium 
with many complicated and non-linear effects, including streaming plasma 
instabilities and electromagnetic field generation. Hence, it is crucial that 
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plasma simulations of internal GRB plasma shocks include lepton-photon 
interactions. 



Solar corona and the solar wind 

Space weather (defined as the interaction of the solar wind on the Earth) is 
in high focus for several reasons. Not only is the Sun our closest star, pro- 
viding us with invaluable data for stellar modelling, but also coronal mass 
ejections from the Sun potentially have impact on our every day life. The 
strong plasma outfiows from the sun can induce large electrical discharges in 
the Earths ionosphere. This may disrupt the complex power grids on Earth, 
causing rolling blakcouts such as the one in Canada and North America in 
1989. Also high-energy particles can be hazardous to astronauts and airline 
passengers. Computer simulations have provided a successful way of obtain- 
ing insight in these complex plasma physical processes. However, in the solar 
coronal and in the solar wind plasma out to distances beyond the earth orbit, 
difficulties arise in finding the right formalism to describe the plasma. Neither 
a coUisionless model based on the Vlasov equation nor an MHD fiuid model 
provides a adequate framework for investigation. The problem has already 
been studied using three d imensional PIC simula t ions but without t aking 
collisions into account (e.g. iBuneman et al. (|l992[ ): lHesse et all (|200lh . 



The corona of compact objects 

The bulk dynamics of accreting compact obj ects have bee n modelled for 
many years using fluid based simulations (e.g. iBalbusI (|2003l ) and references 



therein). Nevertheless, it has been a persistent problem to extract infor- 
mation about the radiating processes. Furthermore in the corona the MHD 
approximation becomes dubious, just as in the solar corona. The environ- 
ment around a compact object is much more energetic than the solar corona, 
and therefore radiative scattering processes play an important role. Pair 
production is also believed to be abundant. Using our new code it would be 
possible to model a small sub box of the corona. The main problem here - as 
in most numerical implementations - is to come up with realistic boundaries 
for the local model. A shearing box approach may be appropriate, but in 
fact we can do even better. 

The size of a stellar mass black hole is around 10^ cm. In a fluid simulation 
we want to model the accretion disk-compact object system out to hundreds 
of radii of the compact object. The normal approach is to use a non-uniform 
mesh. Nonetheless, the Courant criteria, which determines the time step, is 
still limited by the sound crossing time of the compact object. I.e. the time 
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step is limited by the size of the innermost (and smallest) cells in the mesh. 
The very small time step corresponds to those found in a typical particle 
simulation, where the strict time step arises from the need to resolve plasma 
oscillations. Hence data from an MHD simulation could provide temporally 
well-resolved fluxes on the boundaries of the much smaller sub box containing 
the particle simulation. 

In this sense the particle simulation will act as a probe or thermometer 
of the fluid model. The particle model includes the full microphysics in a 
realistic manner and most importantly includes photon transport. Realis- 
tic spectra could be obtained indirectly from the fluid model, testing fluid 
theory against observations. We have already started preliminary work on 
interfacing fluid models with the old PIC code. 

Pre-acceleration in Cosmic Ray acceleration 

Accepting Fermi acceleration as a viable mechanism for accelerating electrons 
and creating the non-thermal cosmic ray spectrum is still left with some big 
unanswered questions. One is that the Fermi mechanism requires injection 
of high-energy electrons while still keeping a large, low-energy population to 
sustain the magnetic turbulence. Hence, a pre-acceleration mechanism needs 
to be explained. 

The shocks in supernova remnants are believed to be cosmic ray accelera- 
tors. However, the Fermi acceleration process in shocks is still not understood 
from first principles but rely on assumptions on the electromagnetic scatter- 
ing mechanism. PIC codes would seem ideal in exploring the mechanism 
from first principles, since they include field generation mechanisms and the 
back- reaction that the high-energy particles have on this scattering agent. In 
Supernova remnants however, the mean free path for Coulomb collisions are 
comparable to the system and particle-particle interactions cannot be fully 
neglected. 

8.2 Implementation 

Implementing any state-of-the-art large-scale numerical code is a big under- 
taking, and can easily end up taking several man years. We estimate that the 
final version of the next generation code will contain more than 50.000 lines 
of code. Starting in February this year, it has taken us three man months 
to implement the current incarnation of the code, which has already grown 
to approximately 10.000 lines. Besides T. Haugb0lle and C. B. Hededal, the 
development is done together with A. Nordlund and J. T. Frederiksen. For- 
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tunately we have a good tradition and expertise for numerical astrophysics 
in Copenhagen and we have been able to port different technical concepts 
and solutions from our suite of fluid codes and to a lesser extent from the 
old PIC code. The aim is to build an extremely scalable code that is able to 
run on thousands of CPUs on modern cluster architectures and utilize MPI 
as the inter node communication protocol. In this chapter we will not go 
further into technical details. Instead we will put emphasis on the important 
concepts and physics and how we have implemented these. 



8.2.1 Concepts 

The two fundamental objects in a particle-in-cell code are the mesh and the 
particles. We have adopted the solver and interpolation routines from the 
old PIC code to solve the Maxwell equations and find fluxes and densities 
on the mesh. The mesh is used to distribute messenger fields - such as the 
electromagnetic fields - and to calculate volume averaged fluxes and densities 
of the particles The latter are used as source terms in the evolution of the 
messenger fields. The particles really r epresent an ensemble of part icles and 



are often referred to as pseudoparticles iBirdsall and LangdonI ()l99lh or large 



particles. A so-called smoothing kernel describes the density distribution of 
a single pseudoparticle on the mesh. In our implementation the volume of a 
particle is comparable to a cell in the mesh. 



Pseudoparticles with variable weights 

The concept of pseudoparticles is introduced since the "real space" particle 
density easily exceeds any number that is computationally reasonable (i.e. 
of the order of a billion particles). The pseudoparticle charge to mass ratio 
is kept the same as the ratio for a single particle. 

In ordinary PIC codes the weight of each pseudoparticle of a given species 
is kept constant throughout the simulation. The benefit is a simple code and 
an unique identity for each particle. The first is a convenience in the practi- 
cal implementation, the second important when understanding the detailed 
dynamics and history of a single particle. 

Notwithstanding possible conveniences, as detailed below in section l5.2.H 
we have decided to improve this concept to a more dynamical implementa- 
tion where each pseudoparticle carries a individual weight. Particles are then 
allowed to merge and split up when a cell contains too many/few particles, or 
when particles are scattered. The concept is sometimes used in smooth parti- 
cle hydrodynamics (SPH), where different techniques have been proposed for 
the splitting and merging of particles. It is both used to adjust the density of 
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individual particles lB(z>rve et al.l (1200 1|) and in the conversion of gas- to star 
particles in galaxy formation models Icovernato et al. (2004). An important 
quality of SPH is its adaptive resolution capabilities. These are important in 
the description of collapsing self-gravitating systems, ranging from core col- 
lapse supernovae to the formation of galaxy clusters, scenarios where matter 
is collapsing many orders of magnitude, and therefore the smoothing length 
or volume of the individual particles is readjusted accordingly. Consequently 
when sphtting particles, or adjusting the weights in an SPH code it is impor- 
tant to match precisely the spatial density distribution of the parent particle 
to the spatial distribution of the child particles. In PIC codes though the 
spatial size or smoothing parameter of an individual particle is determined 
beforehand by the mesh spacing. This is reasonable since we are not inter- 
ested in adaptive resolution but rather a kinetic description of the plasma 
dynamics. Splitting a parent particle with weight Wp into child particles with 
weights U7* is therefor trivial. The requirements of conservation of mass and 
four velocity together with conservation of the density and flux distribution 
in the box, can all be satisfied by setting 



n 

^p = ^< ep = < ipVp = iivi (8.1) 

i=l 



since the smoothing kernel is determined by the mesh spacing, not the mass 
of the individual particle. 

The merging or renormalization of pseudoparticles requires a much more 
thorough analysis. Up to now we have investigated two schemes, one that 
respects conservation of mass, energy and four velocity by merging three 
particles into two at a time, and one where only mass, energy and average 
direction is conserved by merging two particles into one particle. While 
these schemes probably work well for approximately thermal distributions, 
it will easily give rise to a large numerical heating when considering head on 
beam collisions. We believe it can be improved by first selecting a random 
"merger particle" and then find other particles in the local cell, that are 
close to the merger particle in momentum space. A more radical approach 
is to resample the full phase distribution in a cell every time the number 
density becomes above a certain threshold. Nevertheless, it requires testing 
of different extreme situations to find the optimal method to merge particles, 
and it is still a work in progress. 

To obtain the results, that we present in section 18.31 we ran the code 
without merging of the pseudoparticles activated. 
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Scattering processes and splitting of particles 

In Monte Carlo based particle codes the generic way to compute an interac- 
tion is first to calculate the probability for the interaction Ps, then compute 
a random number a. If a < Ps then the full pseudoparticle is scattered oth- 
erwise nothing happens. This probabilistic approach is numerically rather 
efficient and simple to implement, but it can be noisy, especially when very 
few particles are present in a cell. In large particle Monte Carlo codes the 
typical cell contains up to lO'' particles per species per cell (hence "large par- 
ticle"). In our PIC code typical numbers are 10^ — 10^ particles per species 
per cell, since we need many cells to resolve the plasma dynamics. For our 
requirements the probabilistic approach would result in an unacceptable level 
of noise. For example, in a beam experiment the spectra of the first gener- 
ation of scattered particles may come out relatively precise, but the spectra 
of higher generation scattered particles (i.e. particles that are scattered more 
than once) will come out with poor resolution or require an excessive amount 
of particles. Another well known consequence of the probabilistic approach 
is that for a given experiment the precision goes in the best case inversely 
proportional to the square root of the number of particles used in the exper- 
iment. To increase effective spectral resolution we have instead decided to 



Figure 8.1: To implement the scattering of two pseudoparticles we transform to 
the rest frame of the target particle (shown as red/light gray) and computes the 
probability P(n) that a single incident particle (shown as blue/dark gray) during 
a timestep At is scattered on the n target particles. If the incident particle has 
weight m, then k = P[n)m particles will interact and two new pseudoparticles are 
created. 



Laboratory frame 




Target restframe 




take a more direct approach. For simplicity we will here describe the method 
for a two-particle interaction, and disregard all factors converting code units 
to physical units. For example, the weight of a pseudoparticle is proportional 
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to the number of physical particles in the pseudoparticle. Although, these 
prefactors all represent trivial conversions of units, they must be taken into 
account in the real code. 

Consider a single cell containing a single pseudoparticle (red) with weight 
Wt = n and a single pseudoparticle (purple) with weight Wi = m, where 
n > m (see Fig. 18. 1|) . We first select the red particle as the target, since 
n > m, and the purple as the incident particle. We then transform the 
four velocity of the incident particle to the rest frame of the target particle, 
and calculate the total cross section at of the interaction. Conceptually we 
consider the process as a single incident particle approaching a slab of the 
target particle. The number density of target particles in the slab can be 
calculated from the weight wt as pt = Wt/AV, where AV = AxAyAz is 
the volume of a single cell. Given the number density the probability that a 
single incident particle is scattered per unit length is 

During a time step At the incident particle travels Al = ViAt, and the 
probability that a single incident particle is scattered then becomes 



1 - exp [-PiAl] 

WtCTtViAt 



exp 



AV 



The weight of the incident pseudoparticle is Wi = m. Pseudoparticles repre- 
sent an ensemble of particles. Therefore Ps is the fraction of incident particles 
that are scattered on the target. To model the process we create two new 
particles with weight Wmw = WiPs = k. Given the detailed interaction, we 
can calculate the theoretical angular distribution of scattered particles in 
accordance with the differential scattering cross section. Drawing from this 
distribution we find the momentum and energy of the new scattered particles. 
The weights of the target and incident particles are decreased to Wt = n — k 
and Wi = m — k respectively (see Fig. 18. 1|) . 

Our method faithfully represents the actual physics even for small cross 
sections. However, if all the particles are allowed to interact, the number of 
particles in the box will increase at least proportionally to the total number of 
particles squared. This is potentially a computational run away. Normally 
we will have on the order of up to 100 particles per species per cell, but 
to be computationally efficient we only calculate interactions for a subset 
of the particles in a cell. This subset is chosen at random according to 
an arbitrary distribution we are free to select. If the probability that two 



8.2. IMPLEMENTATION 



111 



particles are selected for scattering in a given timestep is Q then the travelling 
length Al simply has to be adjusted as Al/Q. If this arbitrary distribution 
is chosen cleverly, the particles with the largest cross section are actually 
the ones selected most often for scattering, and everything ends up as a 
balanced manner: We only calculate the full cross section and scattering 
as often as needed, and the computational load that is given to a certain 
particle is proportional to the probability of that particle to scatter. We 
rely on the merging of particles as described above to avoid the copious 
production of pseudoparticles. Every time the number of pseudoparticles in 
a given cell crosses a threshold, pseudoparticles are merged and this way the 
computational load per cell is kept within a given range. 

8.2.2 Neutron decay 

Free neutrons not bound in a nucleus will decay with a half-life a little longer 
than ten minutes. The neutron decays into an electron and a proton and an 
electron antineutrino to satisfy lepton number conservation 

n ^ p + e~ + (8.4) 

The rest mass difference of the process (0.78 MeV) goes into kinetic energy 
of the proton, electron and neutrino. Let the neutron lifetime be r in code 
units. If r is comparable to or less than a typical timestep, then practically 
all neutrons decay in one iteration, and it is irrelevant to include them. If 
r is much larger than the total runtime, the neutron can be considered a 
stable particle (unless the neutron density in the box is much larger than 
the proton- or electron density). If instead r ^ a At where a ~ 100, then 
we can select a fraction / of the pseudoparticle neutrons in each cell and let 
them decay. This is done in an analogous manner to the generic scattering 
process described above in section 18.2.11 The weight of the selected neutron 
is decreased with a factor 

(8.5) 

where 7 is the Lorentz boost of the neutron pseudoparticle and / is chosen to 
give reasonable values for the decrease in the weight. At the same time a pair 
of electron and proton pseudoparticles is created with the same weight. The 
generated particles share the excess mass of the process (where the neutrino 
is neglected for now, but could be included in the future). The momenta are 
selected to give an isotropic distribution in the rest frame of the decaying 
neutron. 



exp 



'-fT 
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8.2.3 Compton scattering 

Here we briefly describe a specific physical scattering mechanism, which have 
already been implemented in the code, namely Compton scattering. 

Compton scattering is the relativistic generalization of the classical Thomp- 
son scattering process, where a low energy photon scatters of on a free elec- 
tron. In the rest frame of the electron, the photon changes direction and 
loses energy to the elec tron, which is set in rnotion. The cross section for 
Thompson scattering is Rvbicki and Lightman ( 1979[ ) 



ar = —rl (8.6) 

where tq = e^/(mc^) is called the classical electron radius. The Thomp- 
son scattering approximation is valid as long as the photon energy is much 
lower than the electron rest mass hu -C mgC^ and the scattering can be re- 
garded as elastic. For photon energies comparable to, or larger than, the 
electron rest mass, recoil effects must be taken into account. Measured in 
the electron rest frame we define ei as the photon energy before the scatter- 
ing, €2 as the photon energy after the scattering and 6 the photon scattering 
angle ()8.2|) . By conservation of energy and momentum one can show (e.g. 




Figure 8.2: Schematic view of the Compton scattering process. Impinging on 
the electron, an incoming photon with energy ei is scattered into the angle 9 with 
energy €2- In the initial rest-fi'ame of the electron, the electron will be recoiled to 
conserve energy and momentum. 



Rvbicki and Lightman ( 19791 )) that 



£2 



ei 

1 + -^(1 - cos( 



(8.7) 
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The differential cross s ection as a function of sca t tering angle is given by the 
Klein- Nishina formula Klein and Nishina (|l929> ): lHeitl'e3 (|l954[) 



dac 


A 4 ( ei 




dVL 


2 el \e^ 


e\ 



sin^^ 



(8.8) 



The Klein-Nishina formula takes into account the relative intensity of 
scattered radiation, it incorporates the recoil factor, (or radiation pressure) 
and corrects for relativistic quantum mechanics. The total cross section is 
then 



1 + x f 2x{Wx) 



2x 



- ln(l + 1x) 



where x = hvj {m(?\ 



1 , , , 1 + 3x 

— n 1 + 2x) - 

1x ^ ' 1 + 2x2 



19) 



8.3 Preliminary results 

To test the new code and it capabilities in regard to the inclusion of colli- 
sions, we have implemented and tested a simple scenario involving Compton 
scattering. 




Figure 8.3: 3D scatter plot of a photon beam (black) passing through a cold 
pair plasma (gray). Left panel show initial setup where a photon beam is injected 
in the upward direction. Right panel shows how photons are scattered on the 
electron-positron pairs 
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In the test setup, we place a thin layer of cold electron-positron pair 
plasma in the computational box. From the boundary, we inject a monochro- 
matic beam of photons all travelling perpendicular to the pair-layer (Fig. 18.31 
left panel). As the beam passes through the plasma layer, photons are scat- 
tered (Fig. 18.31 left panel). 

For each scattered photon we sample the weight of the photon and its 
direction (remembering that all particles are pseudoparticles that represent 
whole groups of particles). Fig. 18.41 shows the theoretical cross section as 
function of scattering angle compared with the result from the simulations. 
Four plots for different energies of the incoming photon beam are shown. We 
find excellent agreement between the simulation results and the theoretical 
predictions. 




Figure 8.4: The theoretical Compton scattering differential cross section. We 
have performed a test experiment with an incoming laser beam on a very cold 
electron population. Over plotted the differential distribution is the theoretical 
curve according to Eqs. (|8.8|) and ()8.7|) . 
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8.4 Discussion 

A next generation PIC code that includes many different Icinds of scattering 
process is under development. It will enable us to target problems that 
reside in the grey zone between the MHD and coUisionless plasma domains. 
This domain covers many astrophysical scenarios of great interest counting 
internal shocks in gamma-ray bursts, solar flares and magnetic substorms, 
compact relativistic objects, supernova remnants and many more. 

The concept of splitting/merging particles and keeping individual weights 
of each particle carry many important features. Variable weights represent 
the true statistics of a scattering process in an optimal way compared to the 
Monte Carlo approach. Also, for MPI-parallelization it is crucial that the 
number of particles per cell is kept more or less constant to ensure an optimal 
CPU load-balancing. To localize calculations we are employing a sorting 
algorithm that maintains neighbouring particles on the mesh as neighbours 
in memory. This is not only good for parallelization, but also makes all 
computations very cache efficient; a crucial requirement on modern computer 
architectures. 

To test the infrastructure of the new code we have implemented Compton 
scattering as a simple scattering mechanism. The preliminary results are very 
promising in form of excell ent agreement w i th the theoretical prediction. We 
note that a recent paper bv lModerksi et all ^200^ provide an interesting test 



suite for various kind of particle-photon interactions that can be tested in 
the future. Merging particles has not been satisfactorily implemented yet. 
Parallelization of code is still not there yet, and is necessary to obtain the 
capability of performing truly large-scale experiments. In summary: Work 
has still to be done before we can start to investigate non-trivial astrophysical 
scenarios, nevertheless solid progress has already been made 

This chapter has been written jointly by Christian Hededal and Troels 
Haugb0lle, reflecting the fact that the development process of the next gen- 
eration PIC code has been highly team based. Essentially everybody has 
contributed time and effort to every single source file of the code. It would 
not make sense to write the chapter separately, essentially repeating each 
other and reusing the same figures. 



Chapter 9 



Overall Discussion and 
Conclusions 



The main source of information we have about gamma-ray bursts (GRBs) 
is from multi-wavelength observations of GRB afterglows. The radiation 
from GRB afterglows is generated in collisionless shocks between the external 
plasma (ISM) and the GRB jet. To make the right conclusions about burst 
parameters in these shocks, it is crucial that we have a firm understanding of 
how the radiation is generated. In the current working model, the radiation 
is believed to be emitted from power-law distributed electrons in a magnetic 
field. Following our lack of knowledge about the plasma-physical details, it 
is a general approach to parameterize this ignorance with the dimensionless 
parameters e^, eg and p. As usual, and ee express the fraction of the total 
internal energy that is deposited in magnetic field and electrons and p is the 
slope of the supposedly electron power-law momentum distribution function. 

How the electrons are accelerated and what the origin and nature of the 
magnetic field in the afterglow shock is, remain open questions. 

The main goal behind the work presented in this thesis has been to expand 
our knowledge about the microphysics of collisionless shocks. Attacking the 
problem from first principles, I have used particle-in-cell (PIC) simulations 
to investigate different aspects of these shocks. The results may be catego- 
rized into three groups: Magnetic fields in collisionless shocks, non-thermal 
electron acceleration and radiation from GRB afterglow shocks. I emphasize 
that these three components are highly interconnected and should be seen 
as pieces of the same puzzle. 
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9.1 Magnetic fields in collisionless shocks 

Since the main focus is on ORB afterglows, I have investigated only electron- 
proton plasma shocks. The three-dimensional plasma PIC simulations illus- 
trate from first principles a number of fundamental properties of collisionless 
plasma shocks. 

• In unmagnetized or weakly magnetized relativistic shocks (below mil- 
ligauss), the Weibel two-stream instability will unavoidably create a 
magnetic field that amounts to ~ 10% of the equipartition value (e^ ~ 
0.1) but varies greatly through the shock. 

• The nature of the magnetic field is predominantly transverse to the 
plasma flow, but a parallel magnetic component is also present. The 
parallel magnetic field strength is roughly 10% of the transverse com- 
ponents. An electric field is also present with an energy density that 
amounts to roughly 10% of the magnetic field energy density. 

• The dominating instability is highly non-linear. The electrons are 
rapidly thermalized while the ions form current channels. These cur- 
rent channels are the main source of the magnetic field. 

• Being Debye shielded by the hot electrons, the ion current channels 
are relatively stable and can penetrate deeply into the shock. Earlier 
thinking and concerns were that t he generated mag netic field can only 



survive over an electron skin depth (jGruzinovl 120011 ) . In the simulations 



presented here, the generated magnetic field is sustained for several 
hundred ion skin depths. 

The structure of the magnetic field is highly patched. The transverse 
coalescence scale is comparable to the ion skin depth. A spatial Fourier 
decomposition of the magnetic field shows that the structures follow a 
power-law distribution with negative slope. 

I find that in two-dimensional simulations of the Weibel two-stream in- 
stability, both the evolution of the plasma density profile and the gener- 
ation of electromagnetic field are in agreement with three-dimensional 
simulations. 

With two-dimensional simulations, it has been able to capture a full 
relativistic, collisionless shock. The shock consists of three segments. 1) 
An upstream foreshock with great anisotropy, consisting of instreaming 
ISM plasma and upscattered shock particles. In this region, a strong 
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magnetic field generation is taking place with = 5 — 10%. 2) A 
dense thermalization region where the electrons and ions are close to 
equipartition (eg = 50 — 80%). In this region, the magnetic field is 
relatively weak (e^ ~ 1%) but still strong enough to account for most 
GRB afterglow radiation estimates. 3) A hot downstream region where 
the magnetic field is of the order percents of equipartition. This region 
is clumpy, with large structures that are of the order 1 — 20 ion skin 
depths in size. 

• If the ejecta is already strongly magnetized with a transverse mag- 
netic field (e.g. carried from the progenitor), the resulting dynamics 
is very complex: Particle defiection, charge separation and a mixture 
of the electromagnetic Weibel two-stream instability and the electro- 
static Buneman instability make the resulting field very complicated 
with local reversal of the ambient field. Electrons are accelerated to a 
non-thermal distribution. In the case of strong parallel magnetic field, 
all instabilities are effectively damped. 

• For an ambient magnetic field strength comparable to the ISM- field, 
the Weibel instability grows unhindered. 



9.2 Non-thermal electron acceleration 

There has been a general acceptance of Fermi acceleration as the mechanism 
that produces the desired electron distribution function. Nevertheless, this 
has never been proven in self-consistent numerical simulations. The idea of 
Fermi acceleration in coUisionless shocks is currently facing several problems 
(see Chapter [H for details): 

• It is assumed that the particles scatter on electromagnetic waves, but 
the model does not self-consistently account for the generation of these 
waves nor for the back-reaction that the high-energy particle distribu- 
tion has on the electromagnetic field. 

• Relativistic Fermi acceleration requires the downstream magnetic field 
to be stro ngly turbulent on scales cora parable to the gyro-resonant 
scale fe.g.. lOstrowski and Bednarzlliooi ). This is in contrast with the 



standard synchrotron radiation scenario for GRB afterglows where it 
is assumed that the magnetic field is n early constant on scales much 
larger than the gyro-radius (see however iMedvedev and Chapter 

El). 
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Fermi acceleration requires pre-acceleration to above a certain thresh- 
old. How this pre-acceleration works and how large a fraction of the 
particles that participate in the non-thermal tail remains unexplained. 
In the most well studied mildly relativistic plasma shock (in the Crab 
nebula), obser vations are contradicting wha t is normally assumed about 
this threshold ( Eichler and Waxman . 2005| ). Furthermore, Chandra im- 



ages of the close by Supernova remnant SN 1006 fail to detect an x-ray 
halo in front of the shock, which is to be expected if the standard 
Fermi diffusive shock acceleration theory is correct. One expects an X- 
ray halo around the shock because hig her energy electro ns are expected 
to diffuse further ahead of the shock ( Long et al. . 20031 ). 



I do not claim that Fermi acceleration cannot take place, but urge the need 
for a self-consistent treatment of the problem when computational resources 
allow this. 

I have presented self-consistent, three-dimensional PIC simulations that 
have revealed a new non-thermal electron acceleration mechanism that differs 
substantially from Fermi acceleration. The acceleration is a natural conse- 
quence of the properties of relativistic coUisionless shocks. Acceleration of 
electrons is directly related to the formation of ion current channels that are 
generated in the non-linear stage of the Weibel two-stream instability in the 
shock transition zone. This links particle acceleration closely together with 
magnetic field generation in coUisionless shocks. 

The resulting electron spectrum consists of a thermal component and a 
noi i-thermal com pone nt at high energies. This i s in the line with arguments 
by Rvdd ( 2005[ ) and iRees and MeszarosI ( 2004 ) . In simulations of mildly 



relativistic shocks (P = 3), the non-thermal component vanishes in the ther- 
mal pool. Computer experiment with a bulk flow with P = 15 results in a 
power-law slope p = 2.7 for the electron distribution function. 



9.3 Radiation from GRB afterglow shocks 

I have developed a tool that enables us to obtain radiation spectra directly 
from PIC simulations. The tool has been thoroughly tested and successfully 
reproduces spectra from synchrotron radiation, bremsstrahlung and undula- 
tor radiation from small-angle deflections. The tool has been used to inves- 
tigate the properties of three-dimensional jitter radiation in magnetic fields 
with different turbulent configurations. By tracing ensembles of monoener- 
getic, isotropic distributions of electrons in a random magnetic field (whose 
power spectrum follows a power-law in the Fourier domain Psik) oc k'^) I 
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have computed and studied the resulting spectrum. I have focused on the 
weak field limit (where the defiection angle is less than the beaming cone 
angle, K < 1). For all values of fi, the radiation spectrum has a fiat low- 
energy part (like bremsstrahlung) . For fi < 0, the high-energy part of the 
spectrum follows a power-law with a slope a ^ /i — 1 independent of the 
electron energy. For /i > 0, the fiat spectrum continues to higher frequencies 
than for /i < 0. The spectrum has a hard cut-off for high frequencies. 

For all values of fi and K < 1, the spectrum shifts in frequency with 
7^ and scales with the amplitude of the magnetic field squared B^, in good 
agreement with Larmor's formula for radiated power. In the limit of large 
defiections {K > 1) the spectrum eventually converges through the wiggler 
spectrum to the ordinary synchrotron spectrum. 

I have furthermore examined the radiation from electrons moving in 1) a 
magnetic field generated by the Weibel two-stream magnetic field and 2) a 
magnetic field that has the same average energy density and spectral power 
spectrum as the two-stream generated field. I find that even though the two 
spectra peak at the same frequency, there is a large difference above the peak- 
frequency. The reason is that in the magnetic field from the PIC simulations 
there exists a phase correlations resulting in higher magnetic peak values but 
smaller filling factor. The peak frequency of synchrotron/wiggler radiation 
scales linearly with the maximum magnetic field strength. Furthermore, from 
other PIC simulations we know that the high-energy electrons are found near 
the peaks in magnetic field. 

In simulations of a coUisionless shock that propagates with F = 15 through 
the interstellar medium, the resulting radiation spectrum peaks around lO^^Hz. 
Above this frequency, the spectrum follows a power- law F oc with 
P = 0.7. Below the peak frequency, the spectrum follows a power law F oc u"' 
with a ~ 2/3. This is steeper than the standard synchrotron value of 1/3 
and more compatible with obs ervations. Both t he slope and the peak is con- 



sistent with observations (e.g. |Panaitescijl200ll who finds (3 = 0.67 ± 0.04, a 



peak at 3 X lO^^Hz and F > 10 for the afterglow of GRB 000301c after five 
days) . 

I stress the following interesting point: If one hides all the real phys- 
ical details of the magnetic field in the dimensionless parameters e^, eg, 
and p the conclusion from standard synchrotron radiation theory would be 
that the slope of the electron distribution A^(7) oc is found by solv- 
ing —j3 = —{p — l)/2 = —0.70 — i> p = 2.4 ( which is consiste nt with 
standard analysis of the GRB 000301c afterglow, IPanaitescu 2001). This 
would appear to be consistent with Fermi acceleration and the supposedly 
universal p ~ 2.2 ± 0.2: Case solved! But a closer look reveals that the 
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acceleration mechanism is not Fermi acceleration. The electrons are instan- 
taneously accelerated and decelerated in the highly complicated electric and 
magnetic field near the ion current channels. Strong radiation is produced 
in this process. The electron distribution behind the radiation is a mixture 
of a thermal component for low energies and a power-law component with 
p = 2.7 for high energies ( Hededal et al. . 20041 ) (see also Chapter E}. The 
paths of the electrons are more random rather than circular and the electron 
distribution function varies with shock depth, as does the magnetic topology 
and strength. I repeat and emphasize that the electrons that produced the 
spectrum in Fig. 17.221 are not Fermi accelerated but a natural consequence 
of the Weibel two-stream instability. 



9.4 Future work 

In this thesis I have shown, from first prinicples, that from the collision 
of two plasma shells that are initially cold and non-magnetized, a strong 
small scale magnetic is generated (e^ ~ 0.01 — 0.1), particles are acceler- 
ated to non-thermal distributions and the emitted radiation is quite consis- 
tent with observ ations (comparing Fig. 17.221 and data from GRB 000301c, 



Panaitescijl200l[ ). The next step is to determine the polarization in the syn- 
thetic spectra from the PIC simulations. Even though time did not allow 
this in this thesis, it will not be difficult with the tools that have already 
been developed. With the inclusion of polarization it will be possible to in- 
vestigate how the radiation spectrum and polarization change as functions 
of the viewing angle relative to the jet propagation direction. 

In the near future, computational resources will reach a level where one 
can resolve full three-dimensional collisionless shocks. This will enable us to 
make composite spectra from different simulations with bulk Lorentz factor 
varying with viewing angle, and in turn allow us to test many interesting 
aspects with regard to the jet structure and polarization predictions from 
shock-generated electromagnetic fields. 

Finally, the new generation PIC code that is under development allows 
more of physics to be included in the simulations. Some of the main features 
are: High-energy photons are treated as particles, several scattering pro- 
cesses are included (e.g. Coulomb and Compton scattering), parallelization 
with MPI, pseudo particles with individual weights, etc. With this code it 
will become possible to study also the prompt GRB phase, where Compton 
scattering and pair processes are hkely to be important. 
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Appendix A 



The Mean Free Path in a Blast 
Wave 



In this appendix I make an estimate of the mean free path for Coulomb 
coUisions of a relativistic electron with momentum 7meVe on ions in a plasma 
with density n. 

First we examine the Coulomb collision between an electron and an ion. 
Without loss of generality, we may assume that the electron is travelling in the 
yz-plane along the z-axis and that an ion is positioned in (x, y, z) = (0, —b, 0). 
The ion is surrounded by an electric field. Because of Lorentz contraction and 
symmetry arguments we can assume that the electron will only be affected 
by the component that is transverse to Ve, n amely Ey. In the reference fra me 
of the electron, this component is given by ( Rvbicki and Lightman . 19791 ): 



1 q'jb 



^. = 7^ .,... J (A.1) 



Here t is the time, centered so that the electron is in (0, 0, 0) at t = 0. 
The force felt by the electron is F = qEySy. The change in the electron's 
momentum 6p is found as: 

6p= f Fydt = f qEydt = , ^ (A.2) 

J " J 47reo 6^72t;2t2 + 52 

The pulse from the ion is felt by the electron in the short time interval 
T ~ 6/ (7^,;). Inserting this into Eq. IA.2I we find 

5p=-^ (A.3) 
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APPENDIX A. THE MEAN FREE PATH IN A BLAST WAVE 



We are interested in collisions that alter the impinging electrons momentum 
significantly with 5p ^ •jnieVe, and we can thus find the distance be for such 
a collision: 



6p ~ '-frrieVe 
1 g2 



(A.4) 



47reo \/2^mevl 
Thus, the cross section for the collision is 

ae = Txbl = . (A.5) 

167reQv27^m2?;4 

The collision frequency is Ve = naeVe and from this we find the mean free 
path for a full collision: 

^ Ve nOe nq^ 

In reality, the mean free path somewhat shorter because of accumulation 
of small angle defiections. We can correct f or this by intro ducing a correction 
factor 1/ In A, which is of the order of 0.1 ( Spitzei . 19621 ). 



For a electron in a blast wave that is expanding with Lorentz factor 
7 ^ r = 10, into an interstellar medium with density n ~ lO^m"'^, the 
mean free path for a collision is larger than lO^^m. Comparing this number 
with the typical size of a GRB blast wave ~ lO^^m we conclude that it is 
reasonable to neglect collisions between the ejecta and ISM, and that the 
shock between the two is to be regarded as a collisionless shock. 



Appendix B 
Integration of 



We wish to integrate Eg. 12.81 over all directions n through the solid angle dQ. 



dPr. 



rad 



Hoq^c 



n X Un - /3) X /3 



167r2 



n.f3y 



(B.l) 



Without loss of generalization we define a Cartesian coordinate system with 
the a;-axis along the particles velocity vector f3 and the acceleration vector 
in the x-y plane (Fig. IB.l| . 




Figure B.l: We define a coordinate system with the x-axis along the velocity 
vector j3. For the integration over the unit solid angle dfl around the unit vector 
we use the polar coordinates {(,,^)- The particle is accelerated at some angle 9 to 
the velocity vector. 
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Integration over the unit solid angle dil then yields: 



rad 



rad 



^0 



dn 



sin ^d^dif. 



(B.2) 



With the definition of the polar angles in Fig. IB. II we can writte out /3, /3 
and n: 



(3 = P\ 





n 



cos^ 
sin ^ cos ip 
sin^ sin ip 



:b.3) 



dPrad _ /ipg^C 

dn 



167r2 



(3 sin^ [{(3 — cos^) cosy^sin^ + cos 6^ sin ^] 
+ $ sin ^ (cos 9 cos ^ + cos ip sin ^ sin ^) sin ip 
+ P cos^{{P — cosC,) sin 6' + cos 6^ cos sin 

2 

—$ sin 6* sin^ ^ sin^ (/? 



:B.4) 



Inserting Eq. IB. 41 into Eq. IB. 21 we perform the integration in two steps: 
Jo Jo 



rad 



rad 



dn 



sin ^ dC, dip 



/3V{(11 + 6/3^ + 2/32cos(20 - 2cos(2e) 







(1 + 3(3^ + (3 + f3^) cos(20) - 32/3 cos ^ sin^ i 
sin^ - sin(30}/(16(l -pcos^f) sin^d^ 

2 

f^oQ C 6^2 



127r 



7^/3' (2-/3^(1 -cos(2^)) 
7^/3' (l-/?^in2^). 



(B.5) 



Appendix C 

Spectral distribution of 
Synchrotron Radiation 



Here we derive the spectral distribution of synchrotr on radiation. We follow 
the sta ndard derivation found in many textbooks (e.g. lR vbicki and Lightman 
(197^), but keep the full angular dependency. We start with the expression 
for the radiated energy pr. unit frequency pr. unit solid angle (Eq. 17. 8p 



dujdVt IGTT'^eoc 

Integration by parts and using the relation 

n X [(n - /3) X /3] _ d 
(1-/3 -11)2 ~ Jt 

we can rewrite into the expression 

d^W 



r n X [(n - /3) X /3] ,^(j,_n . v{t')/c),. 
(l-/3.n)2 



n X (n X /3) 
l-/3-n 



dujdVL IGTT'^eoC 
We now expand the terms in the integrand 



/oo 
n X (n X /3)e^'^(*'-^ " 
-oo 



t' - n . r(t')/c) 



i! ^ cos Q sin(/3ct'/ a) 
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/ Br'c^ co-s,Q ,3 
(1 - /5 cos + ——^ i! 



6r| 

, /5^c2cos6' ,3 



(C.l^ 



(C.2) 



(C.3) 



(C.4) 
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Figure C.l: We define a coordinate system with the j;-axis along the velocity 
vector (3. For the integration over the unit solid angle dil. around the unit vector 
we use the polar coordinates {0,Lp). The particle is accelerated at some angle 9 to 
the velocity vector. 



where 6} = 2{1 



[3 cos 6). 
n X (n X /3) = 



j3 [sin 6* cos(/?ct'/ri)e_|_ — sin(/3ct'/r2,)e||] 
P sin 9e± — jS"^ ct' / r Le\\ 



The radiation has to polarization components 



dudVt dudVt dudQ 



(C.5) 



(C.6) 
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duodVL 
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dujdfl 

We now define 
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.(C.8) 



(C.9) 
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dudfl 
dudfl 



exp 



3 1 



dy {C.IO) 



K 



yexp 



3 1 



dy 



(C.ll) 



We now focus solely on the integrals. We hide all the coefficients by intro- 
ducing the variables x and z. 



X = -Y 

2 



a = - cos ep-^x- 



(C.12) 



First, we expand the exponential functions in a geometric representation, 

exp [i{xy + ay^)~\ = cos{xy + ay^) + i sm{xy + ay'^). (C.13) 
The solutions to the integrals of Eq. \C.10\ and Eq. Kllll are expressed by 



Abramowitz and StegunI ()l964f): 



exp [i{xy + ay^)] dy 



yexp [i{xy + ay^)] dy 

'CO 

d f°° 

— / cos(xy + ay^)dy 
dx J_^ 



cos{xy + ay^)dy 



2(3a)-^/VAi [{3a)-^/^x] (C.14) 



(C.15) 



ysm{xy + ay^)dy 
2(3a)-2/3vrAi' [{?>a)-^l^x] 



where A\{z) is the Airy function. The Airy function and its derivative can 
be expressed in terms of the modified Bessel function of the second kind: 



Km 



IT 



-7r^Ai'(z), 

z 



(C.16) 
(C.17) 



where z = (fC)^''^- I our case, z = (3a) ^^^x and by comparing these two 
expressions for z and reinstating the definitions of a and x we find the relation: 



c 



X 



(C.18) 



cos( 



Thus, the Airy function and its derivative in Eq. IC.14I and Eq. IC.15I can be 
replaced with: 



Ai [(3a)-i/=^x] = Ki (xN cos 613 
Ai' [{?,aY^'^x] = -K2 [x/^cosOp^ 



z 

3^ 

z 



(C.19) 
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with z = (3a) ^^'^x. With these expressions we arrive at the final result of 

Ea. irnni and Ea. irnTi 
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q oj 
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